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ABSTRACT 


A dual potential formulation for numerically solving the Navier- Stokes equations 
is developed and presented. The velocity field is decomposed using a scalar and 
vector potential. Vorticity and dilatation are used as the dependent variables in the 
momentum equations. Test cases in two dimensions verify the capability to solve 
flows using approximations from potential flow to full Navier-Stokes simulations. A 
three-dimensional incompressible flow formulation is also described. 

An interesting feature of this approach to solving the Navier-Stokes equations is 
the decomposition of the velocity field into a rotational part (vector potential) and 
an irrotational part (scalar potential). The Helmholtz decomposition theorem allows 
this splitting of the velocity field. This approach has had only limited use since it 
increases the number of dependent variables in the solution. However, it has often 
been used for incompressible flows where the solution scheme is known to be fast 
and accurate. This research extends the usage of this method to fully compressible 
Navier-Stokes simulations by using the dilatation variable along with vorticity. 

A time-accurate, iterative algorithm is used for the uncoupled solution of the 
governing equations. Several levels of flow approximation are available within the 
framework of this method. Potential flow, Euler and full Navier-Stokes solutions are 
possible using the dual potential formulation. Solution efficiency can be enhanced 
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in a straightforward way. For some flows, the vorticity and/or dilatation may be 
negligible in certain regions (e.g., far from a viscous boundary in an external flow). 
It is possible to drop the calculation of these variables then and optimize the solution 
speed. Also, efficient Poisson solvers are available for the potentials. 

The relative merits of non-primitive variables versus primitive variables for solu- 
tion of the Navier-Stokes equations are also discussed. 
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1. INTRODUCTION 


1.1 Background 

The topic of this final report is a particular computational approach for solving 
the Navier-Stokes equations. The Navier-Stokes equations are usually associated 
with the field of fluid mechanics. Solutions to these equations with the appropriate 
boundary conditions model fluid motion. 

An analysis of fluid motion requires the solution for the physical laws of nature: 

1. Conservation of mass 

2. Newton’s second law of motion 

3. Conservation of energy 

These laws can be formulated mathematically, with the help of some assumptions, 
to become the Navier-Stokes equations. Formally, the Navier-Stokes equations refer 
to the mathematical representation of Newton’s second law. It will be more con- 
venient for the purposes here to let the term “Navier-Stokes equations” include the 
representation of all three physical laws above. Assumptions in the development are 
that the coefficients of viscosity are related by a factor of —(2/3) according to Stokes’ 
hypothesis and that the fluid is Newtonian (Schlichting 1979). A Newtonian fluid 
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is one in which the fluid shear stress is linearly proportional to the rate of strain. 
Additional relationships are included as necessary to describe certain processes or 
fluids. Examples are the equation of state for a perfect gas, Sutherland’s formula for 
viscosity and Fourier’s law of heat conduction. 

As one can imagine, the successful solution to the Navier-Stokes equations can 
help immensely in engineering design and optimization. The numerical solution of 
the Navier-Stokes equations can be a complement to experimental and theoretical 
fluid mechanics. Unfortunately, the Navier-Stokes equations are coupled and highly 
non-linear. Only a few exact analytical solutions are available for simple conditions. 
In most configurations of practical interest, numerical techniques must be used to 
obtain a solution. 

Much progress has been made in obtaining numerical solutions to the Navier- 
Stokes equations. Several mathematical formulations for the Navier-Stokes equations 
have been developed. They can be divided into two classifications: 

1. Primitive variable methods 

2. Non-primitive variable methods 

As the name suggests, primitive variable methods solve the Navier-Stokes equa- 
tions using the primary variables as the unknowns. The primary variables are velocity, 
total energy (or a variable related to the energy) and pressure or density. One way 
to think of the primitive variables is that they are physical quantities which one can 
measure in the laboratory. Non-primitive variables, on the other hand, are mathe- 
matically derived variables. They are derived from the primitive variables. The non- 
primitive variables used in this report will replace the primary variable of velocity. 
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The replacements used will be vorticity and dilatation. One can devise techniques 
to measure vorticity and dilatation experimentally, but the direct measurement of 
such quantities is uncommon. Both primitive and non-primitive variable methods 
have been used to obtain solutions to the Navier-Stokes equations by numerical tech- 
niques. Primitive variable methods are the most widely used for three-dimensional 
simulations. Either primitive or non-primitive variables are used for two-dimensional 
flow solutions although most applications of non-primitive variables have been for 
incompressible flows. The following sections in this chapter will discuss some primi- 
tive and non-primitive variable solution methods. The focus of this study will be on 
a particular non-primitive variable method that is extended to compressible viscous 
flow. 

1.2 Primitive Variable Methods 

Numerical methods of solving the Navier-Stokes equations using the physical 
variables have attracted much attention. Several popular techniques will be men- 
tioned here. The solution method depends on whether the flow is incompressible 
or compressible, because the Navier-Stokes equations have a different mathematical 
classification depending on the compressibility. For an unsteady incompressible flow, 
the governing equations are elliptic/parabolic in time. For an unsteady compressible 
flow, the equations are hyperbolic/parabolic in time. 

The most common primitive variable solution method for incompressible flow 
problems involves the use of a Poisson equation for pressure in place of the continuity 
equation. An algorithm which employs this solution method is the SIMPLE (Semi- 
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Implicit Method for Pressure-Linked Equations) procedure (Patankar 1975, 1981). 

Another primitive variable solution method is the artificial compressibility ap- 
proach which modifies the continuity equation to include an unsteady term related to 
pressure (Chorin 1967). The resulting equations are a mixed set of hyperbolic/parabolic 
equations which can be solved using a time-dependent approach. This approach ap- 
plies in two and three dimensions and can be modified to compute unsteady flows. An 
available computer code that uses this method is INS3D (Kwak et al. 1986; Rogers 
et al. 1987; Rogers and Kwak 1988). 

A compressible flow solution is often obtained using a time-dependent or time- 
dependent-like approach. Most schemes utilize implicit methods, for example, the 
Beam and Warming (1978), Briley and McDonald (1977), or MacCormack (1981) 
methods. An available code for these applications is F3D (Steger et al. 1986). Addi- 
tional discussion and references on primitive variable solution methods can be found 
in Holst (1987). 

1.3 Non-primitive Variable Methods 

Methods which in some way replace the velocity with derived variables will be 
discussed here. At the highest level of approximation, potential flows are typically 
! solved using either the velocity potential or the stream function. Examples of their 

use are found in most fluid mechanics textbooks (Currie 1974). By definition, a 
potential flow is irrotational so that the velocity field can be defined by the gradi- 
ent of a scalar function. This scalar function is called the velocity potential. It is 
analogous to the electric field potential. For an incompressible potential flow, the 

i 

i 
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only additional constraint is mass conservation. The equation of mass conservation 
is a Laplace equation for the velocity potential which is easily solved. The velocity 
potential is applicable in two and three dimensions. If the stream function is used 
it is defined to satisfy mass conservation and the Laplacian of the stream function 
satisfies the irrotationality condition. The stream function as defined here only exists 
in two dimensions. In both of the above schemes, the momentum equation (vorticity 
transport equation) is satisfied automatically since the vorticity is zero everywhere 
for the assumption of irrotational flow. 

The vorticity /stream function approach is widely used for solution of the two- 
dimensional incompressible Navier-Stokes equations. This method is also discussed 
in most fluid mechanics textbooks and is treated extensively in the book by Roache 
(1972). 

Many would consider this the limit of practicality for non-primitive variable 
methods. However, there are at least two other noteworthy approaches to solving 
the Navier-Stokes equations in non-primitive variables. Both are valid for two- and 
three-dimensional unsteady flows. These methods are known as the 

1. Vorticity /velocity approach 

2. Vorticity /vector potential approach 

These two schemes will be briefly described and then the focus will be placed on 
the vorticity /vector potential method. The topic of this thesis will cover the vor- 
ticity /vector potential method. This method is also referred to by the the aliases 
scalar/vector potential, vorticity/potential, or dual potential method. Since both 


6 


vorticity and dilatation are used in this work to replace the primitive variable mo- 
mentum equations, it does not seem appropriate to identify the method as the vor- 
ticity/vector potential approach. Instead, the term dual potential will be used here 
following Chaderjian and Steger (1985). This terminology identifies the method as 
one which uses two potential functions in a velocity decomposition. 

1.3.1 Vorticity/velocity approach 

In this method for incompressible flow, the momentum equation is replaced by the 
vorticity transport equation. Derivatives of the vorticity definition then yield Poisson 
equations for the velocity when the continuity equation is used to make appropriate 
substitutions. A more general derivation of the Poisson equations is to take the curl 
of the vorticity and substitute in the vorticity definition from velocity. The identity 
for the vector triple product then yields Poisson equations for the velocity. The 
earliest use of this method was by Fasel (1976). He studied the stability of two- 
dimensional boundary layers using a coupled and iterative algorithm. Dennis et al. 
(1979) used the vorticity/velocity method in the calculation of the cubical driven box 
problem. Orlandi (1987) solved high Reynolds number flows over a backward facing 
step. Other works using the vorticity/velocity formulation are Osswald et al. (1987), 
Guj and Stella (1988), Gatski et al. (1982), Fasel and Booz (1984) and Farouk and 
Fusegi (1985). There have been no reported compressible flow applications of this 
method. However, the dilatation could be used as a dependent variable, as was done 
in this research, to extend the vorticity/velocity method to compressible flow. 
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1.3.2 Dual potential approach 

Derived variables which can be used to represent the three-dimensional continuity 
and momentum equations for an incompressible flow are vorticity, a vector potential, 
and a scalar potential. This is one possible three-dimensional extension of the more 
familiar two-dimensional vorticity /stream function approach. This approach and oth- 
ers that use the vorticity as a dependent variable are appealing because vorticity is 
generally located near boundaries in high Reynolds number flows and subsequently 
diffused and convected away. For a three-dimensional incompressible flow the usual 
procedure in the dual potential method is to solve the vorticity transport equation, a 
vector Poisson equation and a scalar Poisson equation. These equations are derived 
from the continuity and momentum equations where the velocity is defined as the 
curl of a vector potential plus the gradient of a scalar potential. The existence of 
these potentials is easily shown for an incompressible flow since the velocity field is 
divergence free (Aziz and Heliums 1967). 

There has been only one reported formulation of the dual potential method for 
three-dimensional compressible, viscous, unsteady flows (Morino 1985). He derived 
a set of equations for density, vorticity, entropy and the potentials. There have been 
no reported calculations using Morino’s formulation. 

1.3. 2.1 Applications of the dual potential method The dual potential 
method has been applied to inviscid and viscous flow problems. Inviscid flow applica- 
tions include the work of Rao et al. (1989) and Giannakoglou et al. (1988). Rao et al. 
(1987) developed a three-dimensional inviscid rotational flow solver based on the dual 
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potential method. They incorporated a boundary layer interaction scheme for viscous 
flow problems. Giannakoglou et al. (1988) compute two-dimensional steady rotational 
transonic flows in arbitrarily shaped ducts and plane cascades. They decomposed the 
mass flux vector into two potentials. 

In the viscous regime, the dual potential method has been applied to problems 
of three-dimensional natural convection in enclosures (Mallinson and De Vahl Davis 
1973) and three-dimensional incompressible flows in ducts (Wong and Reizes 1984). 
External viscous flows have been computed by Davis et al. (1989). No attempts 
have been reported on the use of this method to solve three-dimensional unsteady 
compressible viscous flows. 

The dual potential method was first applied to natural convection problems by 
Aziz and Heliums (1967). They used the dual potential method to transform the 
Navier-Stokes equations. The transformed equations were solved using an alternating 
direction implicit (ADI) scheme for the parabolic part of the problem (temperature 
and vorticity transport equations) and a successive over-relaxation (SOR) method 
for the elliptic portion (vector potential equations). They tested their technique by 
applying it to the classical problem of convection in fluid layers bounded by solid 
walls in both two and three dimensions. 

Aziz and Heliums showed the dual potential method to be faster and more accu- 
rate than solutions obtained using the primitive variable approach. In fact, though 
the equations are fewer in number for the primitive variable approach, Aziz and 
Heliums report that they are much harder to solve than the equations in the dual 
potential method. The difficulty arises from the highly non-linear nature of the pres- 
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sure equation and the coupling due to pressure in the momentum equations (as in 
incompressible flow problems). 

The technique developed above was used by Ozoe and co-workers (Ozoe et al. 
1976, 1977, 1979, 1985) in solving a variety of natural convection problems. In the 
1985 paper, the problem of three-dimensional turbulent natural convection in a cubi- 
cal enclosure was solved using a two-equation model for turbulence. 

Applications of the dual potential method to incompressible duct flow (through- 
flow) have not been wholly successful due to confusion over the appropriate vector 
potential boundary conditions. The earliest work in this area was by Aregbesola and 
Burley (1977). They presented a numerical finite-difference solution for the equations 
of motion of a steady laminar incompressible flow in two and three dimensions using 
the dual potential method. Wong and Reizes (1984) presented a dual potential for- 
mulation for unsteady incompressible flows in ducts of constant but arbitrary cross 
section. They showed that the method is capable of handling flows over a wide range 
of Reynolds numbers and imply that it can deal with flow situations in which other 
models become inadequate. The dual potential method guarantees a zero divergence 
of velocity while the usual primitive variable method can at best approximate global 
continuity. That formulation was limited to simply connected domains. In a later 
paper Wong and Reizes (1986) showed how to use the dual potential method to solve 
for the three-dimensional flow in multiply connected regions such as annular geome- 
tries. Yang and Camarero (1986) used body fitted coordinates with the dual potential 
method to simulate incompressible laminar flows in a square elbow and in a twisted 
square elbow. The dual potential method in this paper is shown to be applicable to 
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general duct flow situations in simply connected regions. Hafez et al. (1987) used a 
finite element method to solve the steady two-dimensional Navier-Stokes equations 
in a dual potential formulation for subsonic and transonic flows. They computed 
laminar and turbulent flow cases. 

Some viscous external flow solutions were obtained by Rao (1987) and Davis et al. 
(1986) for flow over two- and three-dimensional troughs. Rao (1987) used interacting 
boundary layer theory to supply the vorticity to a dual potential code for the inviscid 
rotational part of the flow. Davis et al. (1986, 1989) use a viscous dual potential 
method for the entire flow field. Matching between the outer inviscid flow and the 
inner viscous region is automatic in their case. 

One possible extension of the dual potential method to three-dimensional com- 
pressible viscous flow has been formulated by Morino (1986). There are no reported 

results in the compressible viscous regime. 

Compressible viscous unsteady flows have been solved by a closely related method, 
however. El-Refaee et al. (1981) used a non-primitive variable method that replaced 
the momentum equation with vorticity and dilatation transport equations. The veloc- 
ity field was obtained from the vorticity and dilatation field by an integral represen- 
tation. They solved compressible unsteady flows and demonstrated that the solution 
field for vorticity and dilatation can easily be limited in their integral representation. 
In this report a similar equation set is used, but the velocity is decomposed into two 
potentials and solved completely by finite differences. The proposed extension of the 
dual potential method to compressible flow would be directly applicable to the vor- 
ticity/velocity method. That is, the dilatation variable would be included to account 
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for compressibility. 

In view of the short list of references on the dual potential method it is evident 
that this method has not been widely implemented in computations. The main 
reasons for this have been the need for the solution of several additional variables as 
compared to the primitive variable approach and the inability to analyze the numerical 
solution process for convergence. There are no good model problems to guide the way. 

Three-dimensional flow solvers are computationally demanding and the introduc- 
tion of additional variables inevitably increases the computer memory requirement. 
However, with the increasing memory of today’s computers, and since Aziz and Hel- 
iums (1967) have shown that the dual potential approach can lead to faster and more 
stable convergence than for primitive variable formulations (for certain problems), 
the vector potential will perhaps play an increasing role in the solution of complex 
three-dimensional fluid dynamics problems (Wong and Reizes 1986). Certainly this 
kind of formulation deserves continued investigation. 

1.3. 2. 2 Advantages and disadvantages Relative advantages of the prim- 
itive variable method and dual potential method are cited in Morino (1985) and 
Richardson and Cornish (1977). The major advantage of working in primitive vari- 
ables is the relative simplicity of the equations and the fact that the primitive variables 
have direct physical meaning. 

The advantages of the dual potential method are: 

1. The vorticity (and dilatation for compressible flow) need only be resolved in 


distinct regions. 
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2. Continuity is automatically satisfied for incompressible flow. 

3. The equations are weakly coupled (at least for inviscid flow). 

4. Good numerical solution routines exist for Poisson equations. 

5. Matching between an inviscid region and viscous region occurs automatically 
because of the velocity decomposition into rotational and irrotational parts. 

Disadvantages of the dual potential method applied to a three-dimensional com- 
pressible unsteady flow are: 

1. The dual potential method involves ten dependent variables whereas the prim- 
itive variable method involves only five to represent conservation of mass, mo- 
mentum and energy. (In two dimensions the number of dependent variables are 
six for the dual potential method and four for the primitive variable method.) 

2. The equations for the dual potential method are more complex than the equa- 
tions associated with the primitive variables (or, they are simply unfamiliar). 

3. The potentials do not have direct physical significance. 

In addition to the natural disadvantages of the dual potential approach listed 
above, there is a lack of available software as compared to primitive variable solution 
methods. The extension of this approach to unsteady compressible viscous problems 
is uncharted territory. 
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1.4 Scope of the Present Study 

It has been the goal of this research to extend the capability of the dual po- 
tential method to compute unsteady compressible viscous flows. An algorithm has 
been developed to provide two-dimensional full Navier-Stokes simulations. A three- 
dimensional algorithm has been developed for incompressible flow only. Test cases 
were computed to verify the ability to compute flow fields ranging from full potential 
flow to flow fields requiring the full Navier-Stokes equations. It has been demonstrated 
in this work that the calculation region can be limited for vorticity and dilatation, 
thus providing a speed advantage for certain flows. 

Several two-dimensional test cases will be presented to test various aspects of the 
dual potential method. Both incompressible and compressible flows will be computed. 

Incompressible flows will be studied for steady, irrotational, inviscid conditions 
and for steady, rotational, viscous conditions. The steady, irrotational, inviscid test 
case is that of flow over a biconvex airfoil (or a bump on a wall). Steady, rotational, 
viscous conditions are simulated for a channel inlet and laminar boundary-layer case. 
Heat transfer calculations will be made for the channel cases with constant wall 
temperature and constant wall heat flux boundary conditions. 

For compressible flow, steady and unsteady, irrotational, inviscid flows will be 
computed and also steady, rotational, viscous flows. The irrotational, inviscid flows 
are biconvex airfoil cases. The steady, rotational, viscous flows are channel inlet 
and boundary-layer cases. The channel inlet flows are computed with constant wall 
temperature and constant wall heat flux boundary conditions at a Mach number of 
0.1. Calculations of the flow over a flat plate are made for a subsonic and supersonic 
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freestream. 

In three dimensions only steady, viscous, incompressible channel inlet solutions 
were obtained. A summary of the test cases to be presented is given below. 

I. Two-Dimensional Cases 

A. Incompressible flow 

1. Steady irrotational inviscid flow 

a) bump cases 

2. Steady viscous flow 

a) channel inlet with and without heat transfer 

b) boundary layer 

B. Compressible flow 

1. Steady irrotational inviscid flow 

a) bump cases 

2. Unsteady irrotational inviscid flow 

a) bump cases 

3. Steady viscous flow 

a) variable property channel flows 

b) boundary layer 

II. Three-Dimensional Cases 

A. Incompressible flow 

1. Steady viscous flow 
a) channel inlet 

B. Compressible flow 

Progress in this research area has not been easy. There is very little guidance 
in the literature on how to proceed with a full Navier-Stokes implementation of a 
non-primitive variable method. The governing equations in non-primitive variable 
form are unfamiliar. Non-linear terms were simply lumped into the source term and 
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the system was solved uncoupled in an iterative manner. As a result, it was necessary 
to employ rather simple test cases to check various aspects of the formulation. 

The work completed here is primarily in the development and evaluation of the 
dual potential method as a flexible approach to solving the Navier-Stokes equations. 
Several features of the method have been highlighted. For example, the solution do- 
main for vorticity and/or dilatation may be limited to certain regions. Also demon- 
strated is the flexibility of the method to accommodate several approximations of 
the full Navier-Stokes equations. This effort has advanced the understanding of the 
dual potential method in the viscous compressible regime. It represents the first 
application of this method to compute throughflow problems with heat transfer. Sev- 
eral basic problems are solved to check out aspects of the algorithm and computer 
code. Only Cartesian grids are used for the test problems. Further evaluation and 
optimization of the method reported herein are left for future work. 

1.5 Organization 

The main body of this report consists of six chapters and two appendices. The 
presentation follows the logical development of the method from equation derivation 
to boundary condition determination, grid generation, numerical algorithm selection 
and, finally, flow simulations. In Chapter 2, the dual potential equations are derived 
from the velocity decomposition and non-primitive variable dependent variables are 
selected to represent the usual primitive variable (or pressure-velocity) form of the 
Navier-Stokes equations. In Chapter 3, the numerical representations of the boundary 
conditions are derived and the numerical algorithms are presented. The Cartesian 
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grid stretching is presented also. A comparison is made of the Poisson equation solvers 
since the Poisson equation solution for the potentials can dominate the computation 
time. In Chapter 4, the solution strategy is outlined and then numerical results are 
reported for two- and three-dimensional test cases. The two-dimensional results repre- 
sent cases from potential flow to situations requiring the full Navier-Stokes equations. 
The three-dimensional results are for incompressible cases only, but are representative 
of the speed of this approach for incompressible problems. Chapter 5 includes the 
overall assessment of this method. Chapter 6 gives some incentives for future work 
on the dual potential method. 

The appendices contain equations for the full three-dimensional Navier-Stokes 
implementation of this method. Also, alternative non-primitive variables are intro- 
duced which could be useful for some problems. 
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2. DUAL POTENTIAL FORMULATION 
2.1 Introduction 

In this chapter the mathematical equations which model fluid flow are presented. 
These equations are the Navier-Stokes equations. They state the conservation of mass, 
momentum and energy for a Newtonian, Stokesian fluid. The usual form of these 
equations has the primitive variables (p, V ,Ei) as the primary unknowns. Using the 
Helmholtz decomposition theorem, the velocity field can be split into a rotational 
part and an irrotational part. Each part is represented by a potential function. 
This decomposition yields a non-primitive variable formulation for the Navier-Stokes 
equations. 

2.2 Governing Fluid Dynamics Equations 

The following equations apply to a continuum fluid. 

The conservation of mass is stated 

Dp -*• -» 

-T + pV • V = 0 

m F 

The conservation of momentum (Newton’s second law), with the assumptions 


( 2 . 1 ) 
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that Stokes’ hypothesis holds and that the fluid is Newtonian, is written 


J)V -» z* d 

P— = F ~ V P + 


Dt 


dx 


J L 


du: duj \ 2 duL 

P \ ^ I - - Si — - 


dxj dxi) 3 tJ dx k \ 


( 2 . 2 ) 


Where F = Xi + Y] + Zk is the body force and is the Kronecker delta function: 


1 if i = j 
0 if i ^ j 

In the energy equation only internal and kinetic energy will be considered im- 
portant. The conservation of energy is then written 

+ P V • V = V -~q + * (2.3) 

where represents heat energy production by external agencies, q is the heat 
conduction and $ is dissipation. Fourier’s law of heat conduction will be assumed so 



9 = 


-kVT 


The dissipation function for a Newtonian fluid in a Cartesian coordinate system be- 
comes 


$ 


2 + Vy + w 2 ) + (v x + Uy) 2 + (t Vy + v z ) 2 + ( u z + to*) 2 

— — + Vy + Wz} ] (2-4) 


The ideal gas equation of state and a viscosity law are used to close the system for 
laminar flow. Constant specific heats are assumed throughout. Reference conditions 
are selected to non-dimensionalize the equations. Reference quantities will be denoted 
by the subscript r. Fluid properties for air will be used in the calculations. 
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For an ideal (or perfect) gas the following relationships exist: 


p = pRT 

e = c v T h = cpT 7 = — R = cp — cy 

c v 

Sutherland’s law of viscosity is used in the form 

j>3/2 


p = C 1 


t + c 2 


where Cj and C 2 are constants for a particular gas. For air at moderate temperatures 
(approximately 200K-1000K), Cj = 1.458 x 10~ 6 [kg/(m s\/K)] and C 2 = 110.4K. 
Power law variations for p were also available in the computer code: 

p ( T \P 

- = ( jr) 0.5 < p < 1.0 (2.8) 

A constant Prandtl number is assumed and thermal conductivity is obtained 
from the definition, Pr = Typical values of the fixed quantities chosen for air 


7 = 

1.4 

’ . • 

(2.9) 

R = 

287 

L k 8 ’ K J 

(2.10) 

Pr = 

r* - o.7 

ky* 

(2.11) 


A reference length will be designated by Lp [m] and a reference velocity by Ur |^j . 
The reference length is taken to be a characteristic length of the problem such as the 
hydraulic diameter for internal flow or chord length for external flow. The reference 
velocity is taken to be the magnitude of either the inlet velocity for internal flows or 
the freestream velocity for external flows. 
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The following reference conditions were set for incompressible flow calculations 
with the Reynolds number and reference length specified: 


= Cl 


t 3/2 

1 r 

Tr + C2 


Pr^p 

Pr 

1.22 

Re/ir 

pr Lr 

U r 

VW 


w 


Lm • °CJ 

kg 


T r = 288.15K (59°F) 

p r 

kr 
pr 

U r 

Mr 

Pr — pr RTr 

m‘ 


kg 


m • s 




m 


N 1 


Reference conditions for compressible flow calculations with the Reynolds num- 
ber, Mach number and reference length specified are as follows: 


T r 

p r 

kr 

U r 

Pr 

Pr 


288.15K (59°F) 


,3/2 


Cl 


Tr + C2 


kg 


m • s 


/ ir c p 
Pr 
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tm^j 
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The following non-dimensional variables are then obtained: 


* 

Lf 


x* = ^- y* = T~ z * = T~ «* = 

Lty* 


z 

Lj* 


(Lr/Ur) 


u = 


u 

~Ur 


* V 

V = 

U r 


w 


- W =TTr » 


p_ 

Pr 


P_ 

Pr 


* 

P = 


PrU? 


T * _ * 


* 

e — 


m 


R* 


R 


(U?/T r ) ° V (i U?/T r ) ^ (U?/T r ) " ^ 

where the variables distinguished by an asterisk are non-dimensional. The non- 
dimensional variables will be used throughout, so the asterisk will be dropped in 
the following. The non-dimensional gas constant above is equivalent to 

1 


c v 


* _ C P 

L T1 


k* = X 


R = 


m2 


so that 




p = pRT 


in dimensional or non-dimensional variables. 


2.3 Derivations 


2.3.1 Velocity decomposition 

The basis of the dual potential method is a splitting of the velocity field into 
rotational and irrotational parts. In this section, the impetus for splitting the velocity 
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field in this way is presented. The Helmholtz decomposition theorem formally permits 
the splitting. 

A useful classification of vector fields is possible using the divergence and curl 
operators (Ames 1977). For V • E = 0 at every point of a region R , the vector field 
~E is said to be solenoidal or divergence free. Physically this means that there are 
no sources or sinks in R. If V x E = 0 at every point in R, the field is said to be 
irrotational. The following classification of vector fields can then be made: 

Class I Solenoidal and irrotational: 

vxi = o v-i = o 

Class II Irrotational but not solenoidal: 

Vxi = 0 v-i^o 

Class III Solenoidal but not irrotational: 

Vxi^o v-i = o 

Class IV Neither solenoidal nor irrotational: 

Vxi/0 V-i^O 


An important theorem in vector field theory is called the Helmholtz decomposi- 
tion theorem. It states that any vector field can be split into a curl free and divergence 
free part. Using the above classifications, it can be observed that the velocity field 
of an incompressible fluid is in Class III and the velocity field of a rotational com- 
pressible flow is in Class IV. Applying the Helmholtz decomposition theorem to the 
velocity vector one obtains: 

V = V <f>+V x~A (2.12) 

It is obvious that the curl free part of the velocity is V <f>. The divergence free 
part (recall from the above that another word for divergence free is solenoidal) is the 
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vector field A . The vector A has been assumed to be solenoidal by design. This is 
an arbitrary but appropriate choice to fit the Helmholtz decomposition theorem and 
to remove the redundancy of describing a three component vector ( V ) using another 
three component vector (A) plus the gradient of a scalar ( <f > ). 

For the classification of the velocity field then, one can compute the divergence 
and curl of Equation 2.12 to obtain: 

V • V = V 2 4> = B (2.13) 

V x V = w = V(V-yl) — V 2 4 (2.14) 

where 4> is the scalar potential, B is the dilatation or rate of volumetric strain, u> 
is the vorticity and A is the vector potential. The vector potential, A , is chosen to 
be divergence free. The Laplacian operator in Equation 2.14 is the vector Laplacian. 
Throughout this report only rectangular coordinates are used so each component of 
the vector Laplacian is similar to a scalar Laplacian. 

The vector potential and vorticity will be represented in three dimensions as 
follows: 


A — A^i -|- A^j 4" A%k 

u> = ui^i + u>2J + u>%k 

The components of the vorticity are obtained from: 


a; = V X V 


(2.15) 

(2.16) 


W 1 = w y ~ v z. 


u >2 = -w x + u z , 


^3 = V X - Uy 
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The velocity components are then: 



\ 

u 


+ -A 3 y - A 2z ^ 

v= 

V 

— 

h + A lz - A 3x 


\ w } 


{ $ z + A 2x ~~ A ly ) 


(2.17) 


In two dimensions only one component of A and u> exist. For the standard 
two-dimensional geometry shown in Figure 3.3 the single components are >13 and 
013. To simplify things in two dimensions, the subscripts will be dropped on the 
vector potential and vorticity so that A and u> refer to the two-dimensional case. The 
velocity components in two dimensions are then: 


v= 


\ 

u 


<f>X + Ay \ 

v ) 


1 

H 


(2.18) 


2 .3.2 New dependent variables 


In the well known incompressible application of the above decomposition, the 

momentum equations become the vorticity transport equation. Continuity is satisfied 

— ► — >■ 

by the solution of a Laplace equation for </>, since V • V = 0 in Equation 2.13. Finally, 
the potentials are used to compute the velocity field. Any other governing equations 
remain unchanged (energy equation, equation of state, etc.). 

For a compressible flow, however, V • V ^ 0 . In this case, V • V = V*<j> = 
B ^ 0 in Equation 2.13. An additional equation is required to give the dilatation, 
B , for the solution of the scalar potential. By counting the number of equations and 
unknowns, one can see that the three-dimensional momentum equations represent 
three equations with three unknowns (or two equations and two unknowns in two 
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dimensions). The vorticity transport equations formed by taking the curl of the three- 
dimensional momentum equations yield just two independent vorticity component 
equations (see Appendix A). In two dimensions, the vorticity is a single scalar quantity 
so again there is one remaining usage of the momentum equations permitted. This 
means that an equation (hopefully for the dilatation, B) can be derived by some 
operation on the momentum equations in either two or three dimensions. Looking at 
possible operations on the momentum equations, one choice is to take the divergence 
of the momentum equation. 

The divergence of the momentum equation yields for possible dependent variables 

either the pressure (as in the pressure Poisson equation) or the divergence of the 

— y — y 

velocity, V • V . The divergence of velocity will be represented by the scalar variable, 
B , known as the dilatation or rate of volumetric strain. 

Another possible combination of the momentum equations gives a scalar variable 
which is the rate of shear deformation (or shear strain rate). Let the shear strain 
rate be represented by the symbol, T. In two dimensions, T = u y - ft;*. This 
dependent variable is formed by taking of the y momentum equation + of 
the x momentum equation. The wall shear stress is simply fiT. The variable set 
of T and u) can form the basis for an interesting computational procedure in two 
dimensions. Unfortunately, one usually hopes to compute the skin friction, not give 
it as a boundary condition. However, this could be a useful inverse type calculation 
procedure (see Appendix B). 
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2.4 Derivation of Non-Primitive Variable Equations 

The dual potential equations will be derived starting from the governing equa- 
tions written in primitive variable form. For purposes of comparison, the number of 
unknowns required in a flow simulation using either primitive variables or the dual 
potentials will be given. The number of unknowns is computed by considering the 
continuity, momentum and energy equations only. 


2.4.1 Two dimensions 

The derivation of the dual potential method in two dimensions will contain the 
fewest simplifying assumptions except, of course, that it only considers two space 
dimensions. The solution capability will be for flows that require the full Navier- 
Stokes equations (i.e., unsteady compressible viscous flow). 

In the primitive variables, the continuity, momentum and energy equations rep- 
resent four equations for four unknowns. These four unknowns may be p,u,v and T. 
These variables are solved using: 


• continuity: 


Pt + (p u h + (P v h I = 0 


( 2 . 19 ) 


• x momentum: 


Du _ dp 1 d 
p ~Dt = X ~fa + itedx 


/ du 2 -* 

1 9 


(du 

dv\ 

K 2 S “ 3 v • y )J 

Re dy 

H | 

[dy 

+ te). 


( 2 . 20 ) 


• y momentum: 

Du =Y _dp 1± 
P T>t 8y Redy 




J_d_ 
+ Re dx 


( du dv 
** \ dy + 8xJ J 


( 2 . 21 ) 
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• energy: 

5I = -( 7 -i)(v.?)T + ^-v.vr + ^ (2.22) 

The pressure may be obtained from the ideal gas law, p = pRT. Any additional 
variables (/t,fc,7,Pr, etc.) must be accompanied by their own equation of state or 
constitutive equation. 

For the non-primitive variables of the dual potential method, the same two- 
dimensional flow requires the solution for and p. The potentials give the 

velocity field. The potentials themselves are determined by the solution of Poisson 
equations derived from the divergence and curl of the velocity (cf. Equations 2.13 and 
2.14). 


To see what is needed for the dual potential formulation, consider the primitive 
variable equations above. They are the governing equations for fluid flow, but now 
it is desired to solve not for velocities directly but rather for the potentials. The 
velocities are subsequently determined from the potentials by solving Equation 2.18. 
The momentum equations must be recast to generate a solution to be used by the 
potentials. Equations 2.13 and 2.14 for the velocity field splitting suggest that the 
divergence and curl of velocity be sought as dependent variables. These are obtainable 
from the * and y momentum equations above by taking the divergence and curl of 
the primitive variable momentum equation. 

The curl of the velocity (vorticity) is obtained as a dependent variable by taking 
the curl of the momentum equation. The group, v* — uy = w, is retained as the 
dependent variable. The divergence of the velocity is obtained as a dependent variable 

Q r\ 

by taking of the x momentum equation and tj— of the y momentum equation and 
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summing. The grouping, Ux + Vy = B, can be retained as the dependent variable. 
The remaining governing equations can be left unchanged from the primitive variable 
equation set. 

The continuity equation can be used to compute the density. The energy equation 
can be used to solve for T or enthalpy, or any other variable that is related to energy. 
The ideal gas equation of state can then be used to compute the pressure if it is 
needed. 

To summarize: The six variables of the dual potential method corresponding to 
the solution of the continuity, momentum and energy equations are determined as 
described below: 

1. The continuity equation is used to compute the density, p. 

2. The curl of the momentum equations gives the vorticity transport equation 
for w. Conservative body forces are eliminated by this operation. Also, for 
incompressible flow, the pressure is eliminated. For compressible flow, however, 
the pressure derivatives remain but can be expressed in terms of other variables 
by using the equation of state. 

3. The divergence of the momentum equations gives the dilatation transport equa- 
tion to be solved for B. 

4. The energy equation is solved for T, or enthalpy, or a related variable. 

5. A Poisson equation is solved for (f> with B as the source term. 

6. A Poisson equation is solved for A with —u> as the source term. 
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The equations for a two-dimensional dual potential formulation will be presented 
below and subsequently solved in non- conservative form. Body forces will be ne- 
glected. 

In two dimensions the dual potential representation of the compressible Navier- 
Stokes equations for constant specific heats are: 



Pt + u Px + vpy 

+ pB = 

0 


(2.23) 

UUJx ~f“ VL Oy 

- (i) 

V 2 u> = 

Si (p,B, 

u,T) 

(2.24) 

Bi -f uBx 

\ + v &y — 

(Si) 

V 2 B = 

S2 (p,B, 

k>,T) 

(2.25) 

Tt + uTx 

+ vT y - j 

U‘„) 

V 2 T = 

S 3 (p,B, 

u,,T) 

(2.26) 




< 

II 

B 


(2.27) 




v 2 a = 



(2.28) 


The source terms for the B and T transport equations are given below. For 
the vorticity transport equation, contains the compressibility and cj contains the 
variable viscosity terms. 

■?! = —Bw + + 3 ^*) “ Px(ux + 3 ^/)] 

+ ~(pxTy - PyT X ) + 


(2.29) 
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1 2 
C 1 = ~ ( u y ^ v x){pxx — Pyy) + ~( v y ~ u x)pxy 


+Px 

+My 


2 1 2 

— (l*>X + By) 7){pxUy + PxVx ~ 2pyUx + —pyB) 

P P 6 

2 1 2 

(— U>y + Bx ) 9(2 pxVy — PyUy — pyVx ~ ~zPxB ) 

P P £ 6 


(2.30) 


For the dilatation transport equation, S2 contains the “compressibility” and c 2 
contains the variable viscosity terms. 


5 2 = ~(«® + Vy + 2ttyV®) - |px(~^y + -B®) + pt/(a>® + -By)j 


BT —2 ^ R 


RT, 2 


— Bv T v p {TxPx + Typy) d 2 “(/>a: + Py) + 


Re 


(2.31) 


c 2 


2 

+ 

+ 


[/xxy(^j/ + i>x) + PxxUx + Myy^y] — ;j~V 2 p 


Px 
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My 
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2 px 


— 2o»y + - 2^-ux - r -Huy + v x ) 

9 3 3 p p p 

2u>® + “ By + Y^B - 2^-Vy - Uy + V®) 


py, 


(2.32) 


The energy equation source term, S 3, contains the compressibility and variable 
thermal conductivity. 


S 3 = -(7 - 1 )BT+ - R ^ — ^B 2 + a > 2 + 4(u®tiy - VyU®)] 


■[fc®7® + ky Ty\ 


(2.33) 


pRePr 1 

The solution strategy for this system of equations is outlined in Section 4.2. 
Briefly, the equations are grouped into an “incompressible” and compressible part. 
An incompressible solution is obtained by computing among the equations in the 
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“incompressible” set of equations given below. A compressible solution requires a 
pass through both the “incompressible” and compressible sets of equations. In that 
case the “incompressible” equations actually contain terms representing compressible 
effects. This grouping is used to allow an incompressible solution to be the starting 
solution for a compressible problem. 


Incompressible < 


vorticity transport equation 

v 2 4> = o 

V 2 A = -w 


Compressible 


dilatation transport equation 
energy equation 
continuity equation 

V 2 <f> = B 


ideal gas law, p = pRT 
property updates : p,k 

The transport equations above are solved using an ADI scheme. Source terms 
of the dependent variable (undifferentiated) are treated implicitly. Derivatives of the 
dependent variable in the source term are treated explicitly so as not to weight the 
off- diagonals. This did provide a slightly faster solution than treating all dependent 
variable source terms implicitly. 

One will immediately notice the many derivatives introduced in the governing 
equations by this method. Even though most of the test cases to be presented are for 
subsonic flows, it was necessary to handle some terms conservatively. In particular, 
the pressure related terms in the dilatation transport equation and the conduction 
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terms in the energy equation are best handled conservatively for the heat transfer 
test cases computed. By treating these terms conservatively it is easier to obtain an 
accurate representation than by expanding out in the chain rule form and differenc- 
ing. 


2.4.2 Three dimensions 

Only the incompressible equations will be derived here. The three-dimensional 
compressible set is given in Appendix A. 

2. 4. 2.1 Incompressible flow The governing equations for viscous incom- 
pressible flow in non-dimensional vector form are: 

• continuity: 

V • V = 0 (2.34) 


• momentum: 


Q Y — ► — y — ¥ — y 1 o — ► ¥ 

+ (V -V) V=- V P+—V Z V +F 
at Ke 


(2.35) 


There are four unknowns to be determined for the primitive variable solution of a 
three-dimensional incompressible flow. Usually these unknowns are u,v,w and p. It 
will be seen that the dual potential method requires the solution for seven variables 
to satisfy the continuity and momentum equations. 

By taking the curl of the momentum equation the pressure is eliminated yielding 
the vorticity transport equation: 

9 w + (y . V) u -{u • V) V = ^-V 2 w + V x F 

He 


dt 


(2.36) 
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If F is a conservative body force (such as gravity) then it too is eliminated by the 
curl. 

Decomposing the velocity vector according to the Helmholtz decomposition the- 
orem as in Equation 2.12 one obtains: 


V = V<f> + V X A 


(2.37) 


Since V ■ (Vx ^4) = 0, the substitution of Equation 2.37 into Equation 2.34 
leads to: 

V 2 <j> = 0 

The relation between the vector potential, A, and vorticity is obtained by taking 
the curl of Equation 2.37. This becomes: 


A - V (V • A) = 


U) 


The vector potential is chosen to be solenoidal so that the above reduces to: 


1 = - « 

Therefore, for viscous incompressible flow, the momentum and mass conservation 
equations and the vorticity definition can be solved using: 


9 " +(F • V) w -(w • V) V= ^-V 2 u + V x ~F 


dt 


Re 


(2.38) 


V> - 0 

V 2 7 = 


(V 


(2.39) 

(2.40) 
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Then, the velocity is decoded from the potentials according to the Helmholtz 
decomposition: 

V =V<f>+ V x A 

If the pressure field is needed it can be obtained by solving for pressure from 
one of the primitive variable momentum equations or by solving the pressure Poisson 
equation. It is also possible to solve for the pressure by computing a force balance on 
an appropriate fluid element since the velocity field and hence the shear stress field 
is already determined by the solution strategy above. 



35 


3. PRELIMINARY ANALYSIS 


3.1 Introduction 

In this chapter the essential parts for assembling a dual potential code are gath- 
ered together. First the boundary conditions for the dependent variables are pre- 
sented. Next, the Cartesian stretched grids are described and, finally, the necessary 
solvers are explained. 


3.2 Boundary Conditions 

A general presentation of the boundary conditions will be given here. The bound- 
ary conditions are applicable for two- or three-dimensional problems. Boundaries have 
been classified as one of the following: 

1. Solid — impermeable boundary (slip or no slip) 

2. Throughflow boundary — boundary crossed by the streamwise velocity 

3. Far-field boundary — a freestream boundary which may be modeled as an 
impermeable boundary, a porous boundary or some other freestream condition 
according to the problem. 
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The boundary conditions are, of course, also identified with a partial differential 
equation. For ease of presentation, only the boundary types as listed above are 
discussed here. To determine what conditions are imposed in a particular problem 
it is necessary to know the classification of the governing partial differential equation 
(hyperbolic, parabolic or elliptic) and the number and type of derivatives of the 
dependent variable. Specific boundary conditions for each model problem will be 
given in figures in the discussion of the results. The solid boundary and throughflow 
boundary conditions for the potentials are thoroughly derived in Wong and Reizes 
(1984) and Hirasaki and Heliums (1970). For the transport variables ( ~u,B,p,T), 
a fully-developed exit condition is specified by dropping second order streamwise 
derivatives (except for density) and upwinding other streamwise derivatives at the 
exit. Example boundary conditions for a two-dimensional channel throughflow case 
are illustrated in Figure 3.3 for reference. 

3.2.1 Scalar potential boundary conditions 

2 

The scalar potential is obtained from the solution of the Poisson equation, = 
B. This is an elliptic equation so a condition on <f> or its derivative must be given on 
all boundaries. 

Since the velocity is decomposed into two potentials, it is useful to ascribe certain 
of the velocity boundary conditions to each potential. It has already been demon- 
strated by Hirasaki and Heliums (1970) that if the scalar potential were used to deal 
with possible throughflow velocities, then simple boundary conditions on the vector 
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potential are possible. Following Hirasaki and Heliums, 

8 4 > . 7t 

— = n .v 

on 

on the boundaries. Thus, the scalar potential has Neumann conditions all around. 
The normal derivative of the scalar potential is the normal outflow velocity at the 
boundary. For a pure Neumann problem such as this, existence and uniqueness of the 
solution are concerns. Existence is ensured if the compatibility condition given by 
Green’s theorem is satisfied. Uniqueness is ensured by prescribing the scalar potential 
at some point. For the test cases to be discussed later, it was possible to make one 
boundary a Dirichlet boundary and still satisfy the Neumann boundary condition by 
virtue of Green’s theorem. The positive n direction in the following is the outward 
normal. 


3. 2. 1.1 Solid boundaries For no flow through the boundary, the condition 

on the scalar potential is 

£ = . 

on 


3. 2. 1.2 Throughflow boundaries An inlet or exit is the best example of a 
throughflow boundary. The condition on <f> is then 

d<t> 


where, 


— = streamwise velocity, say Uj 


“*' = hll 


and Aj is the cross-sectional area of the throughflow boundary. 
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3. 2. 1.3 Far-field boundaries A derivative condition on the scalar potential 
is also used at a far-field boundary as follows: 


— = normal component of velocity through the far-field boundary, usually 0. 
dn 


3.2.2 Vector potential boundary conditions 


The vector potential is obtained from the solution of the Poisson equation, 
y2 A= — u?. This is an elliptic equation so a condition on the components of 
A or its derivatives must be given on all boundaries. 

With the above choice of scalar potential boundary conditions, the boundary 
conditions on A for a simply connected region may be shown to be 


A - 

A ‘= aT = 0 


(3.1) 


where the subscripts t and n denote the tangential and normal components of A 
respectively. An example of the potential boundary conditions for a solid surface in 
the x-z plane are shown in Figure 3.1. 


3.2.3 Vorticity boundary conditions 

Vorticity is only needed for rotational flow computations. It is generated, for 
example, at no-slip boundaries and diffused and convected away. Vorticity may also 
be specified as part of the inlet or initial conditions. Vorticity can also be generated 
by shocks, but such flows will not be considered here. The transport equation for 
vorticity states the conservation of vorticity. Therefore, the boundary conditions are 
extremely important in defining the flow field. 
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Figure 3.1: 3-D Cartesian coordinate system with example boundary conditions for 

the potentials on a solid impermeable surface in the x-z plane 
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x,i index 


Figure 3.2: An initial uniform velocity profile along a viscous boundary for an im- 

pulsive start 

3. 2. 3.1 Solid boundaries The vorticity at a solid boundary is obtained 
from the no-slip conditions. Several different formulations are possible. The wall 
vorticity can be computed from the velocity derivatives or from a Taylor series ex- 
pansion of the vector potential. Consider a wall in the (x, z ) plane. Using the vorticity 
definition from velocities, the wall vorticity (at y = 0) is 


U>1 ='■ Wy 

(3.2) 

u>2 = 0 

(3.3) 

w 3 = ~ u y 

(3.4) 


This method was used by Aziz and Heliums (1967). It was only used in this work 
to compute an initial wall vorticity based on an impulsive start. For example, the 
initial wall vorticity for a two-dimensional impulsive start is given by o> = — uy (see 
Figure 3.2. The finite- difference initial wall vorticity is then: 
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a>(i,l) = - 


u(t,2) 

y(h 2) — y(h i) 


For subsequent global loops, a better approach (Roache 1972) was used to obtain 
the wall vorticity. This better approach uses the vector potential. Consider a vorticity 
producing surface at y = 0 (j = 1) in a two-dimensional flow. The grid may be 
stretched, so let the transverse spacing to the first node point be A; y\ as in Figure 3.2. 
For the i index in the x direction and j index in the y direction, write the Taylor 
series expansion for the vector potential at (i,2) about the wall values as follows: 



dA 

a 1 d 2 A 

/A x2 1 A 


- A ^+Ty 




z,l 


(At/i ) 0 


+ 0(A Vl ) 4 


(3.5) 


From the velocity decomposition in two dimensions (Equation 2.18), u = <f>x + Ay. 
At the wall this is 


/• ,x d< t> 

u(,,1)= ai 


dA 

*, i + d y 


i, 1 

dA 


but, ti(i, 1) is zero for a no-slip boundary. Therefore, 
can be replaced by the scalar potential derivative along the boundary. 


in the Taylor series above 

i,l 


dA 

dy 


i, 1 


dx 


i,l 


(3.6) 


Again using the velocity decomposition in two dimensions, u = <f> x + Ay , the y 
derivative of this equation gives: uy = 4>xy + Ayy. At the solid impermeable wall 
this is: 


du 


_ d 2 <j> 

d 2 A 

dy 

*»1 

dxdy 

*, i + d y 2 
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From the scalar potential boundary condition at a solid boundary, ttH , = 0. Hence, 

u y li, 1 

the second order y derivative of A at the wall is simply: 


&A 

dy 


•M 


du 

u, i _ Q y 

The wall vorticity in two dimensions is the component w 3 = u> — v x — uy. Along the 
wall v = 0, so that u>(i, 1) = — ^ 
the Taylor series for the Ayy term as: 

du\ 


. Therefore, the wall vorticity is introduced into 
i,l 


d 2 A 

dy 2 


i,l 


dy 


i,l 


= -w(i,l) 


(3.7) 


Equation 3.7 is also obtainable from the Poisson equation for the vector potential, 
V 2 A = -u>. At the wall, A = 0, so A X x = 0 and A yy = -w. Substituting for Ay 
and Ayy in the Taylor series of Equation 3.5 one obtains: 


d4> 

A(i,2) = A(i,l)-^ 


i,l 


A 1 . ........ ,0 

Ayi + 2(-^l))(At/ir + 6^3 


hi 


(A n ) 3 + 0(A n ) 4 

(3.8) 


Solving for u>(i, 1) yields the following first order approximation for the wall vorticity 

2 


u>(t,l) = 


(Ayi) s 


A(i, 1) — A(i,2) — 4>x . Ay 1 

hi 


+ 0(Ayi) 


(3.9) 


The vector potential in two dimensions is zero at a solid impermeable wall, so A(i, 1) = 
0 in the above. Using a second order central difference for <f>x with possible stretching 
in the x direction (Section 3.3), the finite- difference formula for the wall vorticity 
becomes: 


u> 


(hl) = 


( AyiY 


u ir ox dx <f>(i + 1,1) -W- l>l) a „ 

A{i,\) - A(i,2) - ^(t,l) g Ay l 


(3.10) 
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In a similar fashion the Taylor series and above substitutions could be carried out 
in the computational plane ( x,y ). The unequal grid spacing is accounted for by the 
metrics to yield the following formula: 


w(t, 1) 




A(i,l)-A(t,2) 


/ >(» + 1> 1) ~ ~ 1>1) 
S!m)V 2 


1 

2 



\ 


1 


(3.11) 


/ 


The above first order form was used most frequently in the results to be reported 
here. A second order formulation also based on the Taylor series is obtainable from 
Equation 3.5 by retaining the ( Ayj) 3 term. The second order boundary conditions 
will be derived here for a three-dimensional case. Consider the x-z plane again. The 
boundary condition on the component will be derived first by starting with the 
Taylor series expansion of the vector potential component Aj. The Taylor series 
expansion for Aj at a mesh point adjacent to the boundary, (y = Ayj at j — 2), is 
given by 


Ai(i, 2, Ac) — Aj(i, 1, k) + 


dA\ 


dy 


i 19^1 

I /» 


a Id 2 A\ 
i,l,k + 


i, l,k 


(Ay x r 


6 dy 3 


i,l,k 


(Ay 1 ) 3 + 0(Ay 1 )^ 


(3.12) 


where the indices i,j,k denote the x,y,z positions, respectively. 

The first term of the expansion is zero since the tangential components of A 
vanish on an impermeable boundary. Hence 


Aj(i, 1, As) = 0 


(3.13) 
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To evaluate the second term, use the velocity component which contains A± y in 
its decomposition. From Equation 2.17 note that 

= 4>z + A 2x - A ly 

Since u» = 0on the boundary, the derivative in the second term of the expansion can 
be written 


dA i 

dy 


i,l,k 


dj_ 

dz 


i, l,k 


+ 


dA 2 

dx 


i,l, k 


(3.14) 


From Equations 2.14 and 3.1, the second order derivative in the third term may easily 
be identified as 

^4r =-w 1 (t',l,fe) (3.15) 

dy L i,l,k 

Here is where the second order method departs from the first order approach given 
earlier. To obtain an expression for the third order derivative of Aj, a linear distribu- 
tion of vorticity and A-^ yy over the first mesh interval is assumed (Wong and Reizes 
1984), so that 

y 


- < * , l(M>k) + “a — [ u, l(*»2>fc) — wi(t,l,fc)] 


d 2 Ai _ d 2 A x v 

r. 9 “ c *> V*> + A. 


d) r - ■ A n 


d‘A 


2 

g y T 


(3.16) 

(3.17) 


Combining the above linear distributions in order to write as a function of A.\ yy 
yields 


(pA i 


Differentiating Equation 3.18 with respect to y yields 


y 


wi(i,2,/c) + -^~2 ~(i,2,k) 


(3.18) 


duj\ d^A\ 1 

dy ~ dy 3 &V\ 


2 

‘"l(>»2, fc) + -—J-(i,2,k) 


(3.19) 
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Solving Equation 3.19 for A \yyy at y = 0 gives 

(3. 

Finally, substituting Equations 3.13, 3.14, 3.15 and 3.20 into Equation 3.12 gives 
a second order relationship between the boundary vorticity and the potentials, that 
is, 

<i(m ,*) = - t. ■ *> + A 2 x (i, i, *)] 

+ ±A lyy (i,2,k)+0( A n ) 2 (3.21) 

Similarly, the z component of vorticity at this boundary is found to be: 

1. 1) = + ~^ A 2z^ 1.*) + \ A 3yy(h 2 . *) + O ( Ay,) 2 (3.22) 

The third and final vorticity component at this boundary is computed from the 
no-slip condition and the definition of vorticity to be 

u> 2 (i,l,fc) = 0 (3.23) 

as already stated at the beginning of this discussion. In actual use, the second order 
vorticity boundary condition increased the CPU time with no noticeable improve- 
ment in accuracy. Roache (1972) also reports that second order vorticity boundary 
conditions can be less stable and less capable than the first order boundary conditions. 



d * A \ ( : i ».\ w l(*,l, fc )“ w l(*> 2 ,*) , 

3 — a i 


dy 


At/i 


Ayi 


d 2 A i 

dy< ' 


|o>i(t,2,fc)-|- ^A-^(i,2,fc)| 


3. 2. 3. 2 Throughflow boundaries The only throughflow boundaries used 
here are the inlet plane and exit plane for two- and three-dimensional problems. The 
inlet flow field may be specified as either rotational or irrotational. This is controlled 
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by the inlet plane transverse velocity. For two-dimensional flow, the inlet v velocity 
is controlled to provide rotational or irrotational flow as follows: 


inlet condition < 


rotational 

irrotational 


; set v = 0 

; compute v from the decomposition 


( 3 . 24 ) 


The inlet plane vorticity is then computed from the definition, lo = V x V . 

For a two-dimensional rotational inlet condition, the inlet transverse velocity is 
set to zero and used throughout the code. The irrotational inlet condition is achieved 
by computing the inlet plane transverse velocity from the decomposition and using 
that non-zero inlet transverse velocity throughout the code. The streamwise velocity 
component is unchanged since it fixes the inlet mass flux. 

At the exit the second order streamwise derivative is dropped and other stream- 
wise derivatives are upwinded. 


3. 2. 3. 3 Far-field boundaries The far-field vorticity must be specified or 

determined from the specified upstream flow conditions. It is not safe to extrapolate 
the vorticity to a far-field boundary, or any boundary since that does not account for 
the physics of vorticity production, convection and diffusion. Far enough away from 
a surface generating vorticity, the vorticity should be zero. If the velocity is known 
at the far-field, vorticity may be obtained from the definition. 


3.2.4 Dilatation boundary conditions 


The dilatation, B , is only needed for a compressible flow. Boundary conditions 
are developed from the expected velocity boundary conditions and the definition of 
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the dilatation, 

B = V • V = u x + Vy -f Wz 
or by applying the continuity equation, 

Pt + u Px + v Py + ™pz + pB = 0 

The treatment of this variable is the most uncertain. There is little guidance 
in the literature. El-Refaee (1981) used the dilatation in his non-primitive variable 
formulation. He extrapolated to obtain the dilatation at the boundary. The boundary 
values of dilatation were relaxed and set to zero as the solution approached steady 
state. El-Refaee solved external flow problems only. 

3. 2. 4.1 Solid boundaries For a viscous flow along an impermeable wall, all 

velocities are zero. This reduces the continuity equation to 

Pt + P B = 0 

For steady boundary conditions on density, or at steady state it is clear that 

B = 0 

Otherwise, the density time derivative or the velocity derivatives must be evaluated 
to compute 

, from continuity (3.25) 

wall 
or, 

B = ux + vy + w z , from the definition of B. (3.26) 

In any case, B = 0 at the boundary at steady state. 
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3. 2. 4. 2 Throughflow boundaries The only throughflow boundaries used 
here are the inlet plane and exit plane for two-dimensional compressible flow problems. 
The inlet flow may be specified as either rotational or irrotational as stated in the 
section on vorticity boundary conditions. From the definition of dilatation, B = 
u x +vy, an inlet Dirichlet condition can be computed. For the rotational flow inlet 
condition (which corresponds to a uniform parallel inlet flow), v = 0 at the inlet. 
This gives B = u x for the inlet boundary condition. For the irrotational flow inlet 
condition, v ^ 0 and in general Vy ^ 0. This gives B = u x -\-vy for the inlet boundary 
condition. It has been found that an inlet boundary condition of B — 0 may be used 
for either inlet condition. By setting B = 0 at the inlet one avoids the ambiguity of 
specifying the dilatation inlet condition using velocities which in turn depend on B. 
At the exit, the second order stream wise derivative is dropped and other streamwise 
derivatives are upwinded. 

3. 2. 4. 3 Far-field boundaries Generally, the velocity field is unchanging at 
the far field, so B is given by the unchanging velocity derivatives. A typical far-field 
condition is uniform parallel flow so that B — 0. 

3.2.5 Density boundary conditions 

The density boundary conditions are set using one of the following approaches, 


depending on the problem: 
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1. 13y applying the boundary-layer assumption for viscous compressible flow at a 
solid boundary and the ideal gas law, 

Py = 0 


so, 


p(h 1) = 


p(t,2)r(i,2) 

T(i,l) 


2. By using the Bernoulli equation for inviscid irrotational compressible flows, 
P= [l + ^M^ (i_„2_ t) 2_ w 2^ ; FT 


3. By application of the continuity equation written along the boundary, 


Pt 4" u px + vpy -f pB = 0 


3.2.0 Temperature boundary conditions 

The temperature boundary conditions are either set as Dirichlet conditions or 
derivative conditions based on a prescribed wall heat flux. First through fourth order 
polynomial derivative conditions are included in the computer code as options. The 
inlet temperature field is user specified. At the exit plane, second order streamwise 
derivatives are dropped and other streamwise derivatives are upwinded. 


3.2.7 Velocity boundary conditions 

The velocity is not a primary variable in this method. Velocity boundary condi- 
tions are dictated by the flow physics and are used to develop boundary conditions on 
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the derived variables, such as the dilatation boundary conditions above. The condi- 
tions to use are straightforward — e.g., no-slip at solid boundaries for a viscous fluid. 
It is important to point out that the velocity boundary conditions are used, but not 
necessarily enforced by the dual potential formulation. For example, a no-slip condi- 
tion is used to develop boundary conditions for the vorticity at a solid viscous wall. 
However, there is a small slip velocity computed by the velocity decomposition. The 
slip velocity goes to zero as the grid is refined. To see how the slip velocity arises, con- 
sider the velocity decomposition for a two-dimensional flow over a flat plate oriented 
as in Figure 3.2. The velocity decomposition for this case is given by Equation 2.18, 
repeated here for convenience. 

\ , , , 

U <px T Ay 

V ) <j>y — Ax 

At the flat plate surface, the boundary conditions on the potentials are: 

4>y = 0 

A = 0 




It is obvious that the v component of velocity will be zero both analytically and 
numerically (provided the same difference formula is used to compute v as was used 
to enforce the boundary condition on <f>). However, the u component will only be zero 
if - - Ay. This is not exactly satisfied numerically. 

3.3 Cartesian Grid Clustering 

Simple independent variable transformations are used to allow for stretching of 
the two-dimensional Cartesian grids (Anderson et al. 1984). The stretching trans- 

I 

I 
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y,j index 

////////////////////////////////////^^^ 

inlet: 
u> = 0, or 

U) = Vx — Uy 

B = 0, or 
B = u x + Vy 
T = 1 
P= 1 
>1 = 0 

<f>X = u(l,j) 

—£ x 

u> = f(A,<f>) ’ 

B = 0 

walls: T = T w or ^ given 

p = fixed by py = 0 and T distribution 

>1 = 0 
<f>y = 0 


exit: 

(j^xx = 0 

Bxx = 0 
Txx — 0 
px upwinded 

A x — 0 

<f>X = u 


Figure 3.3: 2-D channel boundary conditions 
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formations may be readily applied to a third coordinate direction for use in a three- 
dimensional Cartesian grid. The transformations used here are from the family of 
general stretching transformations proposed by Roberts (1971). 

The coordinate transformation maps the physical plane ( x,y ) into the computa- 
tional plane (x,y): 

(x,y) 0 x,y ) 

where, x — f(x) only and y — f(y) only. Also, the transformations are scaled so that 
Ax = Ay = 1. The computational grid coordinate values then correspond to the 
grid indices like (x, y) = (i — 1, j — 1). This simplifies the coding and avoids divisions 
in the numerical algorithm for a slight speedup advantage. In the following, NI is the 
largest x index and NJ is the largest y index in the domain. NI and NJ correspond 
to the maximum dimensions of the physical grid L,h (see Figure 3.4). 

Applying these transformations to the governing fluid dynamic equations requires 
the following partial derivatives: (These have been simplified since x and y are only 
functions of the respective coordinate directions.) 

d dx d 

dx dx dx 

d _ dy d 
dy dy dy 

d 2 _ (dx\ 2 d^_ (dh\ d 

dx ^ \® x ) dx '2 + dx 

cP_ = d_ 

dy 2 \dy) dy 2 dy 
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d ^ dx dy d 2 

dxdy dx dy dxdy 

The above metric terms are computed numerically using second order accurate 
differences for the unequal mesh spacing. Analytical calculation of the metrics are 
possible for the stretchings to be presented. It was found, however, that the numerical 
calculation provided better flow results. 


3.3.1 Clustering near boundaries 


Take the x direction to be streamwise. The transverse direction is then y with 
“walls” at y = 0 and y = h. The following transformation allows packing near the 
inlet and near one or both walls. 


x 

y 


Nin - ln {[°~ + 1 - (*/£)] / k - 1 + (f/x)jh 

1 ln[(<r + l)/(<7-l)] j 

a + 


1 < <t < oo (3.28) 


a _ a > OP + fe (2« + 1) /h] - 2a} / {/3 - [y (2a + 1) /h) + 2a}). _ 

In [(/?+!)/(/?- 1)] 


This is designed so that 0 < z < (NI - 1) and 0 < y < (NJ - 1) for 0 < x < L, 
0 5: V < ^ with Ax — Aij[ = 1. Equation 3.29 for y packs near y = h for a = 0 
and near both walls equally for a = 0.5. The inverse of equations 3.28 and 3.29 are 
needed to construct the physical grid (x,t/). The inverse for the above transformation 
is readily found as: 


j: ( g + l)-^-l)[g±l]( 1 -Nfa) 


x 


(3.30) 
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y = hi - K + Sgn(/c) 


nJ=T a 

{(3 + 2a) ^ ** — 0 + 2a' 

* ,1*1-“, 
!+l|S 1_ “ ) 


(3.31) 


(2o + 1) ( 


where, 


sgn(/c) 


1 for k > 0 


( — 1 for k < 0 

The inverse for y has been modified using k to direct the y clustering to either 
wall as described below. The stretching parameters in the above transformation have 
the following affects: 


x direction clustering 

1 < a < oo — The stretching parameter a clusters more points near x = 0 as a 
approaches 1. The grid becomes more uniform as x — ► oo. An 
essentially uniform x grid is generated for <r = 10. 


y direction clustering 

1 < /3 < oo — The stretching parameter (3 controls the y direction clustering 

(spacing ratio). Maximum clustering is achieved as (3 approaches 1. 
An essentially uniform y grid is generated for (3 = 10. 


The stretching parameter a is either 0 or 0.5. 


a = 0 


k > 0 clusters near y = h only. 
k < 0 clusters near y = 0 only. 


a = 0.5 Cluster points near y — 0 and y = h equally, k makes no difference 

when a = 0.5. 


An example of the above clustering is shown in Figure 3.4. 
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Figure 3.4: Typical 2-D channel grid, (a) Physical plane using <7 = 1.05, a = 

0.5, /3 = 1.2; (b) Computational plane 



This clustering technique was used for the bump problem (isolated airfoil). It 
can also be used to cluster points about an obstacle located within the grid. The 
equations given will work in either the x or y directions and are designed to cluster 
near a single point or to symmetrically cluster about an object by reflecting the 
generated grid about the line of symmetry. 

The x coordinate scheme for the bump problem will be given here. This requires 
an odd number of i points and assumes that the bump or airfoil is always placed in 
the center of the x grid. (The bump in the test cases is simulated by the blowing 
condition rather than occupying x,y space.) 

As stated previously, x is simply given by the grid point index, i - 1, so that 
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(a) (b) 

Figure 3.5: Typical 2-D grid for an isolated airfoil at y = h. (a) Physical plane using 

r = 5.0, a = 0,/c = 0,(3 = 1.5; (b) Computational plane 


Ax = 1. The inverse transformation is: 


x 


x c \ 


sinh 


1 + 


r vni/2 


w 


sinh (tW) 


(3.32) 


where 


W = -In 



, 1 + (^ T - 1 ) {ift). 


0 < r < oo 


(3.33) 


The stretching parameter, r varies from zero (no stretching) to large values which 
produce the most refinement near x = Xc ■ An example of this stretching is shown in 
Figure 3.5 where the grid has been refined near the line xc and reflected about the 
line of symmetry at x = Lj 2. 
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3.3.3 Arbitrary user defined clustering 

The dual potential computer program will accept any user defined clustering 
so long as x = f(x) only and y = f(y) only. These may be input as arbitrary 
x,y points or as equations relating x,y to x,y as in the above examples. The code 
automatically scales the metrics to form a two-dimensional computational space with 

Ax = Ay = 1. 


3.4 Poisson Equation Solvers 

Three different methods have been used to solve the Poisson equations for the 
potentials. One method is the vectorized point Gauss-Seidel scheme with successive 
over-relaxation (SOR). The other two are alternating direction implicit (ADI) type 
schemes (Mitchell and Griffiths 1980). The two ADI type schemes will be distin- 
guished as follows: 

1. A scheme formed by factoring and then splitting into a two step formula similar 
to the D’Yakonov (1963) approach. This will be referred to as the approximately 
factored (AF) scheme in this discussion. 

2. A Peaceman-Rachford type scheme with a half-time step level (Peaceman and 
Rachford 1955). This scheme will be referred to hereafter as the ADI scheme. 

In summary, the three schemes to be used to solve the Poisson equations are: 

1. Vectorized point Gauss-Seidel with SOR 


2. AF scheme 
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3. ADI scheme 

The model equation = S will be used to demonstrate the above algorithms. 

For the AF and ADI schemes, the Poisson equation is written with a fictitious time 
derivative: + V^3> = S. The iterations represent “time” levels with time step h. 

The time step may also be thought of as a relaxation parameter. At convergence the 
time derivative term is negligible. 

For a Cartesian grid with possible stretching according to the coordinate trans- 
formations in Section 3.3, the equation to discretize is: 



The computational plane grid spacing is unity so Ax and Ay do not appear in 
the algorithms below. 


3.4.1 Vectorized point Gauss-Seidel with SOR 


The finite-difference form of the equation to solve is: 

d 2 y \ % 

dy 2 I 2 


9 *) 2 g 2 , ( 8 M % , ( ? V ) 2 s 2 
Tz) 2 +UJ S V + 


^most recent _ g 


Since Gauss-Seidel is a point iterative method, the exact application of the above 
algorithm will depend on the mesh point ordering. The Gauss-Seidel method is based 
on immediate use of the most recent values. Therefore, the < $ mos * recent use( j j n 
the above equation is either or Solving the above equation iteratively 

by points will not vectorize due to data dependency. This can be illustrated by a 
simple example. Consider the five-point formula finite-difference scheme for the two- 
dimensional Laplace equation, V^U = 0, on a uniform Cartesian grid with Dirichlet 
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boundary conditions: 

U i+l,j ~ 2U i ,j + U i-l,j U i,j+1 ~ 2U i,j + U i,j-1 A 

(A *) 2 + (Atf* 

As a further simplification take A® = Ay = 1 so that solving this by points one 
would code the following: 

DO 40 J=2,NJ-1 

DO 40 1=2 ,NI-1 

U(I,J) = 0.25 * (U(I-1,J) + U(I+1, J) + 0(1, J-l) + U(I,J+1)) 

40 CONTINUE 

The dependency of U(I, J) on U(I-1 , J) inhibits vectorization. Notice, however, 
that for a fixed J , the even I indexed values of U on the left hand side depend only 
on the odd I indexed ones on the right hand side. A vectorization strategy is now 
apparent. The data dependency is removed from the computation by “coloring” the 
grid as a checkerboard and updating the U in two sweeps. At a starting J, the odd 
I index points may be thought of as black squares on a checkerboard and the even 
I index points may be thought of as red squares. At the next J, the odd I points 
are then red squares and the even I points are black squares. This red-black coloring 
continues in J until the grid is patterned like a checkerboard. In one sweep the black 
points are updated using only red points and in another sweep the red points are 
updated using only black points. This is easily implemented by incrementing the I 
loop by 2. Some initial work is required to determine the starting and ending I indices 
for each J. Therefore, in two sweeps the solution is iteratively updated and the code 
will vectorize. The compiler, however, will not recognize that the data dependency 
has been removed. The programmer must direct the compiler to vectorize the loops. 
This strategy is coded below. Note that the starting and stopping I indices are a 
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function of J. It is a simple alternating function, such as 2, 3, 2, 3,. . .for IBSTRT(J) and 


then 3, 2, 3, 2, ...for IRSTRT(J). 

C. . .BLACK POINTS 

DO 40 J=2 ,NJ-1 

C... COMPILER DIRECTIVE: IGNORE VECTOR DEPENDENCIES IN THE I LOOP 
DO 40 I=IBSTRT(J) ,IBEND(J) ,2 

U(I,J) = 0.25 * (U(I-1,J) + U(I+1,J) + U(I.J-l) + U(I,J+1)) 

40 CONTINUE 

C. . .RED POINTS 

DO 41 J=2 ,NJ-1 

C. . .COMPILER DIRECTIVE: IGNORE VECTOR DEPENDENCIES IN THE I LOOP 
DO 41 I=IRSTRT(J) , IREND(J) ,2 

U(I,J) = 0.25 * (U(I-1,J) + U(I+1,J) + U(I.J-l) + U(I,J+1)) 

41 CONTINUE 


This may be combined with SOR for a further speed advantage. The exten- 
sion of this vectorization concept to three dimensions is straightforward. The three- 
dimensional problem may be solved as a stack of two-dimensional checkerboards, as 
a four-color point method, or by extending the idea of colored points to colored lines 
and solving by lines rather than points. The red-black strategy and other vectorizable 
structures are discussed in Gentzsch and Neves (1987). 


3.4.2 AF scheme 


Using first order temporal differencing and second order spatial differencing on 
the model Poisson equation one obtains: 


<j>fe+l _ 


+ 


( d2 A *s.+(?£\ 2 &+ 

Tz) 2 + U) y + 


' <Py\ s y 

dy A 


$k+l _ g 


(3.34) 
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It is convenient to use the “delta” form of the dependent variable so define 

To form a “delta” in the difference equation above subtract the spatial derivatives of 
from each side. Multiplying through by the time step (or relaxation parameter), 
h, gives 



The “delta” form allows easy implicit handling of the boundary conditions. 
Steady Dirichlet conditions are automatically satisfied by the fact that A fe $ = 0 
on the boundary. Steady Neumann conditions are easily handled by reflecting the 
A^$ at the boundary. The derivative function cancels so there is nothing to be added 
to the right hand side. The actual Dirichlet values or derivatives are input by the 
source term spatial derivatives of 

An approximately factored form of the equation above is: 



where the factors can be denoted as Xj and L 2 so that the AF representation is 
L x L 2 A fc $ = RHS. 

As it is written, the algorithm is implemented by sweeping implicitly in the x 
direction then in the y direction. The solution (A^$) is attained in the two steps: 
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Step 1: Li TEMP = RHS 

Step 2: I 2 $ = TEMP 


3.4.3 ADI scheme 


The starting form of the difference equation is the same as for the AF scheme 

(Equation 3.34). The ADI scheme splits the calculation for 4> into two steps. In the 

first step, the x derivatives are treated implicitly at the iteration level k + while 

jt 

the y derivative terms are lagged. The provisional solution is denoted as $ 2 . 

The second step obtains the solution, by discretizing the y derivative terms 

implicitly at the iteration level k + 1 and using the x derivative terms at the k + ^ 


level. 





d 2 x 
dx 2 


6 X 


. k-\- 


1 


hJ 2 







hj 




#*+1 = - 

2 




2\dx) 2 


dx 2 ) 2 J * J 


1 

2 


3.5 Poisson Solver Comparisons 


An assessment of the two-dimensional Poisson solvers was made on test prob- 
lems for the scalar and vector potentials. The test problems are from incompressible 
channel throughflow cases. The three solvers were tested on stretched and uniform 
grids with various aspect ratios. The 1-2 norm of the error is used to com P are the 
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convergence history for each method: 

h = \i 


1 _ 2 
number of points 


However, the requirement used for a converged solution is 

|$k+l _ 


dj^+1 

*max 


< e 


( 3 . 35 ) 


where $ represents either potential and k is the iteration level. The tolerance used for 
the comparisons below is e = 10 . This rather strict tolerance can yield a decrease 

in the X-2 error by as much as 6 orders of magnitude. It was found that the scalar 
potential had to be obtained to this accuracy for reliable flow solutions. 

All of the Poisson solvers were coded with the boundary conditions incorporated 
in the solution algorithm. The relaxation parameter used for the vectorized SOR 
method will be designated by u>. The range for u> is 0 < u < 2. The relaxation 
parameter used for the AF and ADI schemes will be designated by h , the fictitious 
time step. The time derivative term which was added to the Poisson equation sim- 
ulates a parabolic problem, $ t + V 2 $ = S. The way this is written, marching is 
only permitted in the negative “time” direction. Therefore, a negative h was used 
to march the AF and ADI solvers. The need for a negative h is also evident in the 
numerical representation of these schemes since a negative h will add to the diago- 
nal term of the coefficient matrix. The AF and ADI schemes were solved using a 
vectorized tridiagonal solver. 


64 


4> = o 

4>y — 0 x 

Figure 3.6: Scalar potential test problem 

3.5.1 Scalar potential test problem 

The scalar potential is given by the Poisson equation, = B. For an in- 
compressible flow the dilatation, B, is zero. The test problem for incompressible 
flow through a channel is then, V 2 </> = 0, with boundary conditions as shown in 
Figure 3.6. The exact solution for this case is <f> = * (notice that the problem is 
actually one-dimensional). This is a difficult problem to solve numerically due to the 
many Neumann boundary conditions. For incompressible problems, whether steady 
or unsteady, the equation for the scalar potential can be solved once and for all. An 
efficient solver may not seem important for the incompressible case. For a compress- 
ible flow, however, the dilatation field will be computed by a time marching method 
so the scalar potential will have to be solved as often as every global iteration. The 
cases to be studied in this report are mostly subsonic and the dilatation may be ex- 
pected to be small. Therefore, the incompressible solution may be used as a starting 
solution for a compressible problem. Also, in some cases the dilatation may be so 
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small as to be treated as a perturbation of the solution to the incompressible prob- 
lem. For these reasons and the fact that the exact solution is easily obtained for the 
incompressible case, the incompressible problem will be used to assess the solvers for 
the scalar potential. 


The first comparison on the scalar potential test problem is with uniform grids. 
The aspect ratio, defined as a = is varied for each case. The results for a 21 x 21 
grid are shown in Table 3.1 and for a 41 x 41 grid in Table 3.2. The SOR method 
has trouble with this problem and requires an unusually high optimum relaxation 
parameter, a>. It can be seen from these results that the SOR method cannot compete 
with ADI or AF when the grid aspect ratio is far from 1. Stretched grids will be 
necessary in the solution of viscous problems so this immediately excludes the SOR 
method for use in solving for the scalar potential, especially when B is non-zero. 
Notice that the AF and ADI schemes solve the problem in the same number of 
iterations regardless of the cell aspect ratio. Also, the optimum relaxation parameter 
for AF and ADI can be reasonably predicted from the results shown in Tables 3.1 
and 3.2 for 0 < a < 1. This is the most likely range of a for channel type viscous 
flow geometries. The optimum h for the AF scheme behaves like 

ha ~ ha* (1 + 64 log — ) 
x a 


The optimum h for the ADI scheme behaves like 



where h a j denotes the optimum h for a uniform grid with a = 1. This can be used to 
get a reasonable estimate for the optimum relaxation parameter to use for a stretched 
grid. 
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Table 3.1: Poisson solver comparison on a 21 x 21 uniform grid 


Uniform 21 X 21 grid 



~ o — 

iterations to convergence 

CPU time (s 

) 


method 



0.01 


HHBU 

0.01 

MFLOP 




7253 

0.0634 

0.349 

4.11 

39.7 

IESH 

I^Kl 


3 

0.0016451 

0.001626 

0.001667 

54.1 

IEK 

138 


138 

0.1234 

0.1265 

0.1236 

64.6 


| Optimum relaxation parameters (a> or h ) 

method 

ESQ 


0.01 

SOR 

1.898 

1.9854 

1.9991 

COM 


-10240 

-1310720 

ADI 

-0.0089 

-0.89 

-89 


Testing these solvers for the scalar potential on a realistic grid with stretch- 
ing gives the convergence behavior shown in Figure 3.7. The grid is stretched as 
shown in Figure 3.4. The cell aspect ratios for this stretched grid range from a = 
0.4044-0.01265. Notice that the AF scheme converges over a very wide range of the 
relaxation parameter, h. In contrast, the SOR method has a very limited range of 
relaxation parameter which gives convergence in the thousands of iterations at the 
very fastest. The extremely good performance of the AF scheme on this test problem 
is incidental. It is explained by the fact that the exact solution for this test case 
is <f> = x, and the factorization error introduced in the AF scheme contains cross 
derivative terms. The cross derivatives and hence the factorization error therefore go 
to zero quickly. The optimum conditions for the solvers on this problem are shown 
in Table 3.3. The convergence history at the optimum conditions for the solvers is 

shown in Figure 3.8. 
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Table 3.2: Poisson solver comparison on a 41 x 41 uniform grid 


Uniform 41 x 41 grid 

mmm 

iterations to convergence 

CPU time (s) 



a =1.0 

0.1 

0.01 

1.0 

0.1 

0.01 

MFLOP 

SOR 

195 

1260 

8105 

0.306 

2.01 

12.7 

55.3 

AF 

3 

3 

3 

0.00468 

0.00475 

0.00483 

72.8 

ADI 

258 

258 

258 

0.719 

0.713 

0.701 

82.1 


Optimum relaxation parameters (u> or h) || 


a =1.0 

0.1 

0.01 

SOR 

1.9482 

1.9925 

1.99932 

AF 

-160 

-10240 

-1310720 

ADI 

-0.0045 

-0.45 

-45 


With stretching, AF and ADI outperform the vectorized Gauss-Seidel easily on 
this problem. With non-zero dilatation and more complicated geometry the AF 
scheme is not expected to display such a tremendous advantage over ADI as in the 
example here. Also note the high relaxation parameter required for the fastest con- 
vergence by the SOR method. 


Table 3.3: Poisson solver comparison on the scalar potential test problem 


41 x 41 stretched channel grid 

method: 

optimum 
relaxation param. 

iterations to 
convergence 

CPU time (s) 

MFLOP 

SOR 

1.99644 

1954 

3.33 

53.7 

AF 

-80000 

3 

0.00509 

70.9 

ADI 

-3.0 

142 

0.4059 

80.8 
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Table 3.4: Poisson solver comparison on the vector potential test problem 


41 x 41 stretched channel grid 

method: 

optimum 
relaxation par am. 

iterations to 
convergence 

CPU time (s) 

MFLOP 

SOR 

1.721 ' 

39 

0.0629 

54.1 

AF 

-0.003 

35 

0.0531 

70.1 

ADI 

-0.00071 

105 

0.2143 

82.2 


3.5.2 Vector potential test problem 

The vector potential test problem has been taken from an incompressible channel 
throughflow case. The equation and boundary conditions are shown in Figure 3.9. 
For these conditions, the SOR method did not exhibit extreme sensitivity to the grid 
aspect ratio and compared favorably with the AF and ADI methods. The Poisson 
equation for the vector potential, with three of the four boundaries having Dirichlet 
conditions, is easily handled by SOR using optimum relaxation parameters in the 
approximate range 1.2 < < 1.8. This equation must be solved as often as every 

global iteration for a rotational flow. 

The vector potential equation was solved here on a 41 x 41 grid stretched as 
shown in Figure 3.4. A converged incompressible channel flow at Re = 300 provided 
the vorticity source term. The behavior of the solvers over a range of relaxation 
parameters is given in Figure 3.10. The optimum conditions are given in Table 3.4 
and the L2 error for the optimum convergence of each solver is shown in Figure 3.11. 
In this case the AF and SOR schemes compare favorably in CPU time. In an actual 
flow solution, however, the ADI scheme has achieved convergence faster than AF or 
SOR when the vector potential solver is called at each global iteration. 
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Figure 3.9: Vector potential test problem 

3.5.3 Summary of Poisson solver experience 

This was not meant to be an exhaustive study of methods for solving Poisson’s 
equation numerically. Two very specific Poisson equation problems were studied 
to assess the numerical solvers. The results for the scalar potential test problem 
are particularly perplexing. The point SOR method performs poorly while the AF 
scheme is unbelievably fast. This is at least partly due to the one-dimensional nature 
of the scalar potential test problem. Point SOR cannot take advantage of the one- 
dimensional nature of the solution because it solves the domain pointwise. In the first 
iteration, point SOR creates a two-dimensional distribution. This happens because 
the solver cannot sense the information from opposite boundaries simultaneously. 
This is unlike the AF and ADI schemes which can immediately detect the influence 
of opposing boundaries in the implicit direction. The AF and ADI performance on the 
scalar potential test problem could be dependent on the sweep order, especially since 
the solution is one-dimensional. This was not tested, however. The biggest advantage 
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Figure 3.10: Convergence behavior of the Poisson solvers on the vector potential test problem 
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Convergence History 
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Iteration number 

Figure 3.11: Convergence history for the vector potential test problem 
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for the AF scheme is that the factorization error is fortuitously zero because of the 
one-dimensional solution to the scalar potential test problem. 

The vector potential problem doesn’t offer the anomalies that the scalar potential 
problem does. Point SOR and the AF scheme are very competitive. The poorer 
performance of the ADI scheme is perhaps due to the fact that it is not written in 
delta form. 

The apparent overall best solver for the Poisson equations in the dual potential 
method is the AF scheme. At the optimum convergence condition it converges in 
the fewest iterations and in the least CPU time. The AF scheme converges rapidly 
over a large range of relaxation parameter. This is important in getting started with 
a solution for a new problem. Almost any reasonable relaxation parameter, h, will 
work whereas some effort is required to get a fast solution by the SOR method. The 
ADI scheme is a consistent performer for the very different problems tested, but its 
optimum is not as fast as AF at its optimum for the conditions tested above. However, 
in actual flow calculations where the solver may be called at each global iteration the 
ADI scheme is always competitive. SOR is competitive when the boundary conditions 
are Dirichlet. It has a narrow band of relaxation parameter giving fast convergence 
as compared to either AF or ADI. This narrow band and steep slope (see Figures 3.7 
and 3.10) is undesirable for predicting a good relaxation parameter to use on a new 
problem. One rule of thumb, however, is that the optimum relaxation parameter for 
SOR increases as the number of grid points increase on the same problem. 

The above results and discussion suggest that the scalar potential Poisson equa- 
tion is best handled by AF or ADI. The vector potential problem is best handled 
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by AF or SOR, though the ADI scheme was found to be competitive on actual flow 
problems. In actual use, the AF scheme was the easiest to code, followed by the ADI 
scheme and then the vectorized point SOR. All three schemes were fastest when the 
boundary condition information was incorporated into the algorithm. Not indicated 
in the simple test cases above is the fact that SOR loses any advantage on large grids. 
The computational effort for SOR is greater than for the AF or ADI schemes. Beyond 
some grid size then, the SOR scheme should be dropped in favor of the more efficient 
AF or ADI schemes. Pointwise iterative schemes tend to bias the solution depending 
on the sweep direction. The effect of sweep direction was not studied for the SOR 
scheme above. 

For the three-dimensional (incompressible) problems the vector potential equa- 
tion was solved using point SOR. The exact solution of the scalar potential was used 
so that the Poisson equation was not solved numerically. 


3.6 Time Marching ADI Solver for the Transport Equations 


The ADI scheme presented here is based on the one proposed by Douglas (1962). 
The transport equations for a *,B,p and T in three dimensions may be written as 


8S_ _ ( s (PS 
dt \ Re dx 2 U dx 


ds\ Is d 2 s ds\ 

+ I — — 2 - t>— I + 


3 d 2 S as . ,, . 

lS^r- w aT l+c5 -‘' 


Re dy 2 dy 

where S represents any of the dependent variables ,B,p or T and represents 
the coefficient of the diffusion terms. The term 0 includes all remaining terms. The 
source term cS may be lumped into 9 , but here it will be treated implicitly. 


The above is written in a uniform Cartesian grid. A stretching transformation can 
easily be introduced, but would only clutter the development here. The stretching 
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transformations were used in Section 3.4 on the Poisson equations to provide an 
example of their use. Here, the complete three-dimensional algorithm on a uniform 
grid will be sufficiently complicated without metric terms included. 

In the development of the ADI scheme for the transport equations, the following 
definitions will be used 


n 


t 

8x 


i Ax, yj = j Ay, z fc = fcAz 
n At 


~$x , /i \ 


u + |u| V x U - jttj Ax 


A* 




Oy 

' Sz = “"li; + (1 ~ o) 


V -f |v| Vy V — |v| A 


2 Ax 

y 


Ay 


Ay 


w + |w| V z w — |tu| A z 


A; 


Az 


(3.36) 

(3.37) 

(3.38) 

(3.39) 

(3.40) 


where 8x,8y and 8 Z are hybrid finite-difference expressions for the convective terms 
an d respectively. The weighting parameter, a, is in the range 0 < 

d < 1. A central difference formula is obtained for a = 1 and an upwind formula 
for a = 0. A different weighting parameter may be used for each direction. Note 
that this formula is written for the physical grid with uniform spacing. If there is 
stretching, the metrics are required. 

The notation in Equations 3.38-3.40 could be confusing. It is important to note 
that the V and A in the numerator are operators and the A in the denominator 
represents increments. 

Let Q denote the finite-difference approximation to S and define: 


Q n = Q n {i,j,k) = Q(x { ,yj,z k ,t n ) 


(3.41) 
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^ (3 - 42) 

^ < 3 - 43 > 

= fe8^)« n < 3 - 44 ) 

An alternating direction method for the solution of S may now be developed. 

From a known solution at time level n, a first estimate of the solution for time level 
n + 1 may be obtained by evaluating half of the differences with respect to x at a 
provisional time level denoted by n + 1* (Aziz and Heliums, 1967). The source term 
cS will also be evaluated at the implicit time level, so that 

^(Q n+1 %Q n )+cQ n+1 * + 6yQ n + 6 z Q n 

= (q u+1 * - /At + 9 (3.45) 

where 9 can be evaluated at n,n + j or n 4- 1. In the algorithm described here, 9 
was evaluated at the n time level only. In the following, one asterisk denotes the first 
approximation with more asterisks for successive estimates. Each successive estimate 
is made by evaluating half of the differences in one direction implicitly. 

^*(Q n+1 *+C? n ) + +<? n )+cQ n+1 ** +6 z Q n 

= [(g n+1 ** -Q n )/At +0 (3.46) 

\s x (c? n+1 * + Q n ) + \s y (Q n+1 ** + Q n ) + l -S z (Q n+1 + Q n ) + cQ n + l 

= [(<? n+1 -Q n )/Af] +9 (3.47) 

The equations 3.46 and 3.47 may be simplified by subtracting Equation 3.45 
from 3.46 and Equation 3.46 from 3.47, respectively. The resulting algorithm, which 



is second order accurate in time and up to second order accurate in space, is: 


( 8 x + 2c — 2/ At) Q n ^~^ — — (8 x + 2 8y + 2 8 Z + 2 / A t'j Q n + 29 ( 3 . 48 ) 

(8 y + 2c-2/At)Q n+1 ** = 8 y Q n + (2c-2/At)Q n+1 * ( 3 . 49 ) 

(S z + 2c-2/At)Q n+1 = 8 z Q n + (2c-2/At)Q n+1 ** ( 3 . 50 ) 

This is the algorithm used for the three-dimensional solution. A tridiagonal 
system of linear algebraic equations for Q is solved three times to obtain Q n "^. The 
first pass treats the x direction implicitly, the second treats y implicitly and the final 
pass treats z implicitly. The alternating direction implicit method used here reduces 
to the one proposed by Douglas if 8 x ,Sy and S z are replaced by 6£,8y and 6%. Also, 
the non-linearities in the method of Douglas are included as part of 0 only. 

In two dimensions the provisional time level is n + The time increment is ^ 
and the entire differences in one direction are taken at the implicit time level. In the 
(x,y) plane, the algorithm to solve the transport equations can be written 

= 8 x Q n+ ? + 8 y Q n + cQ n+ ? -9 ( 3 . 51 ) 

= 6 x Q n+ ? + 8 y Q n+1 + cQ n+1 - 9 ( 3 . 52 ) 

In a form ready to code, the two-dimensional algorithm becomes: 

(2/ A t-8 x -c)Q n+ ? - (8 y + 2/At)Q n -9 ( 3 . 53 ) 

(2/ A t-8 y -c)Q n+1 = {8 x + 2/At)Q n+ ? -9 ( 3 . 54 ) 


- Q n 
At/2 

QU+1 _ + ^ 

At/2 
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4. NUMERICAL RESULTS 


4.1 Introduction 

In this chapter the test cases will be discussed. It will be seen that the dual 
potential method can be used to solve the full potential equation, the Euler equations 
and the Navier-Stokes equations. The results generated from the dual potential code 
will be designated as DP in the comparisons. 

4.2 Solution Strategy 

The solution strategy has been broken down into three parts: an incompressible 
segment, a compressible segment and an iterative update. Certain steps are common 
to both the incompressible and compressible procedure. The incompressible segment 
and iterative update are always needed. The transport equations are solved uncoupled 
and iterated only once at each time level. The ADI scheme for the transport variables 
(u> ,B,T,p) provides up to second order accuracy in time and space. The Poisson 
equations may be iterated to convergence at each time level for a time accurate result. 
It is possible to limit the iterations on the Poisson equations to speed up the solution 
for a steady state case. The iterations may also be limited for a time accurate solution, 
but numerical experimentation is required. 
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The step by step calculation procedure is: 

“Incompressible” segment: 

1. Using known values (or assumed initial values at startup) for p,^o ,B,T, V at all 

points at a given time level, n, solve the vorticity transport equation by an ADI 
method to give u> 

2. Solve the vector potential equation using the n+1 level vorticity. If the calculation 

is for a constant property, incompressible flow without heat transfer, go to step 
7. Otherwise, continue with step number 3. 

Compressible segment: 

3. Solve the dilatation transport equation for B n using n level source terms. 

4. Solve the energy equation for T n +1 using n level source terms. 

5. Solve the continuity equation for p Tl ^ ^ using n level source terms. 

6. Update the properties p, k and p using re + 1 level quantities. 

Iterative update: 

7. Compute the scalar potential, <f>. 

8. Update the velocity to time level re + 1. 

9. Solve for the vorticity boundary conditions at time level re + 1. 


10. Check for steady state convergence and stop if the solution has converged. Oth- 
erwise, transfer the just computed re + 1 level results to the re level and repeat 
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the sequence from step 1 until steady state convergence or to some point of 
interest. 


Though the segments above may be somewhat misnamed (since the equations 
listed under the “incompressible” segment must be solved for a compressible solution), 
the idea is to suggest possible solution strategies. One strategy which may be used 
to obtain a compressible flow solution is to start with a converged incompressible 
solution. The governing equations of the dual potential method are neatly segregated 
by incompressible/compressible or constant property /variable property criteria. 

The convergence requirement is (usually) applied to all dependent variables. For 
example, a compressible case requires that c o , B,T,p, A and <j> all satisfy the conver- 
gence requirement. The requirement is 

_ $Tl| 


* n+1 

^max 


< e 


(4.1) 


applied at each point, where $ represents any dependent variable and n is the itera- 
tion level or time level according to the nature of the variable. The tolerance is usually 
e = 10~ 5 . This type of convergence requirement is more stringent than an L 2 norm 
and more helpful for a self-interrogating scheme to determine where to concentrate 
computational effort and where to avoid computation. In particular, for the variables 
used in this method, the vorticity and dilatation can be negligibly small. For exam- 
ple, in external boundary-layer flows the vorticity and dilatation can be expected to 
approach zero at some distance away from a disturbance. This point by point error 
checking, normalized on a representative maximum field value, can indicate that the 
largest errors are near the disturbance and initiate a check to see if function values, 
such as B or u> , are small at some distance away. If this is the case, the code can 
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automatically reduce the computational field for that variable. The more common £2 
measure of error that other investigators use has been found to be somewhat deceiv- 
ing. Convergence in the £3 sense ma y n °t satisfy convergence in the absolute sense 
above. Since the £3 measurement averages out the error contribution of each point, 
it doesn’t identify particular locations of large error. That is, the £2 measurement 
may satisfy its convergence criterion and yet there may be isolated areas of large and 
unacceptable error. Hence £2 convergence may have little to do with the accuracy 
attained by the solution. 

Regarding the actual tolerance used in calculations, it has been observed that 
the vector potential error tolerance may be loosened to speed the solution with no 
negative effects. The scalar potential, on the other hand, requires a tight tolerance 
in all cases. The potentials, given by solutions of the Poisson equations, can be very 
time consuming to calculate. 


4.3 Two-dimensional Cases 

The dual potential formulation is rather thoroughly tested in two dimensions. 
Example problems will be computed for incompressible and compressible flows. A 
potential flow solution will be computed for the flow about a thin biconvex airfoil (or 
bump on an inviscid wall). An unsteady calculation for the compressible flow about a 
thickening airfoil will also be discussed. Steady viscous flows will be computed for the 
channel inlet case and flat plate boundary layer. The results with discussion follow. 
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Figure 4.1: Boundary conditions for laminar incompressible irrotational flow over a 

bump 


4.3.1 Incompressible flow 


4. 3. 1.1 Steady irrotational inviscid flow 


Bump cases This test case and the corresponding compressible test 
cases for the flow over thin biconvex airfoils were included to demonstrate the calcu- 
lation of simple irrotational inviscid flows. By considering only very thin airfoils, the 
small disturbance theory boundary conditions may be used. The geometry is therefore 
simple so that the dual potential formulation can be tested with fewer complications. 
Also, velocity potential and stream function solutions are readily obtainable for com- 
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parison. The biconvex airfoil may also be considered to be a thin bump on an inviscid 
wall. Airfoil and bump will be used interchangeably in referring to the bump test 
cases. 

Problems of this type were solved extensively in the 1970s as models for the 
flow over helicopter blades (Ballhaus and Steger 1975; Beam and Ballhaus 1975). 
The bump is generated by blowing according to the small disturbance theory (Ashley 
and Landahl 1965). The airfoil in this case is very thin so that the freestream is only 
slightly perturbed by the presence of the body. Tangency is then imposed by requiring 
the resultant flow velocity at the body to be tangent to the thin body. Since the body 
is assumed to be very thin, the tangency condition can be applied at the airfoil chord 
line. This permits solutions to this problem to be made on a simple Cartesian grid. 
The assumptions for this case are steady, irrotational, inviscid, isentropic flow. 

A sketch of the problem is shown in Figure 4.1. This problem may be solved in 
many ways using potential methods. Only one of the potentials in the dual potential 
method is required to solve this flow. Either one of the potentials can be used with 
little modification to the code. Small perturbation approximations will be used in 
this problem. The velocity components will be represented by 

u = U oo 4" u (4*2) 

v — v (4*3) 

where u and v are perturbation quantities. These are all non-dimensional quanti- 
ties normalized on the freestream velocity so that Uoo = 1.0 and the perturbation 
quantities are assumed to be much less than 1. 

The scalar potential, </>, may be used without modification for a potential flow 
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solution to this problem. The definition of the scalar potential is designed to satisfy 
the irrotationality condition automatically: 

<j> = ux + vy (4-4) 

Conservation of mass is then satisfied by the solution of 

V 2 <f> = 0 (4.5) 

The boundary conditions for a solution of this case using the scalar potential are 
shown in Figure 4.1. 

If the vector potential, A , is used it is actually treated as the stream function for 
irrotational incompressible flow. This is a redefinition of the vector potential that is 
used elsewhere in this development. The similarity of the stream function with the 
vector potential is evident in the Poisson type of equation and Dirichlet boundary 
conditions. This permits the single two-dimensional dual potential code to solve 
this problem using either potential. Let A be the stream function in the following 
development. The stream function is designed to satisfy continuity automatically: 

A x = u (4.6) 

A y = -v (4.7) 

The irrotationality condition is then written 

V 2 A = 0 (4.8) 

The boundary conditions on the stream function are all Dirichlet with A = y. For a 
compressible flow solution the vector potential requires a boundary condition modi- 
fication to emulate the stream function. 
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The scalar potential is perhaps the most natural potential function to use for 
this case since there is no need to change its definition. It simply assumes the role of 
the traditional velocity potential. 

Two different bump geometries were computed by the dual potential code. So- 
lutions have been computed for the parabolic arc and the sine wave arc bumps. The 
parabolic arc has stagnation points at the leading and trailing edges. The sine wave 
bump was included to avoid stagnation points in the early development of the com- 
pressible code. The bumps are described using the small disturbance theory boundary 
conditions. That is, the bumps are generated by blowing through the boundary. The 
parabolic and sine wave arcs are described by; 


parabolic arc: 

y 

= 2 TX (1 — X ) 

(4.9) 

sine wave arc: 

y 

PP 

= [1 — cos (27 t;e)] 

u 

(4.10) 


for 0 < x < 1. r/2 is the height of the parabolic arc and PP (for peak-to-peak) is the 
height of the sine wave arc. The sine wave amplitude is but the sine wave arc 
is used so that the bump meets the boundary with zero slope. The bump geometries 
are indicated in Figure 4.2. The computational domain was chosen to be three chord 
lengths in the streamwise direction and two chord lengths away from the airfoil to the 
far field. The airfoil is centered in the streamwise direction. A uniform Cartesian grid 
was used for the results presented here with 61 streamwise and 41 transverse points. 

The airfoil surface will be defined by the flow tangency condition imposed along 
the mean surface of the body. Flow tangency at the airfoil surface can be represented 
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flow over a parabolic arc airfoil 






Figure 4.4: Pressure coefficient for Mqo = 0*6 flow over a sine wave arc airfoil 
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by 


dy 

dx 


surface ^°° “ 


U ( 


oo 


i 1 + Tt) 


(4 - n) 


In the above equation, jp— is small compared to 1 so it may be dropped leaving the 
following approximation for the tangency boundary condition: 


dy_ 

dx 


r>sj 

surface ^°° 


surface 


= <f>y 


for Uoo = 1. 


surface 


(4.12) 


Finally, the tangency condition is applied at the mean surface of the airfoil so that 


dy_ 

dx 


v 


surface 


u< 


oo 


4>y 


for U, 


grid boundary 


oo 


1. 


(4.13) 


Isurface 

This small perturbation type of boundary condition does not require the use of a 
body-oriented grid and is therefore very easy to implement. The above development 
has been for the scalar potential. As used in the test problem, the boundary condition 
on the scalar potential is 

4^ o < x < i 

<py = •( X surface (4-14) 


0 


otherwise 


where the solution domain is -1 < x < 2. 

A similar development for the stream function gives the following small pertur- 
bation boundary condition 


A=l 


y + local airfoil height 

y 


0 < x < 1 
otherwise 


(4.15) 


The airfoil height referred to here is the half thickness of the airfoil at a particular x 
position. 
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The pressure coefficient is defined as 

C P = Tfw < 416 > 

oo ^oo 

where V<£, = Uqo + The pressure coefficient for the bump cases is computed by 
an approximation to the definition 

C P = f^20 K (4.17) 

\pooV& U °o 

This approximation is obtained by dropping squares of the perturbation velocity com- 
ponents. In terms of the scalar potential or stream function, the pressure coefficient 
is approximated by 


C p ss2(l-fe)=2(l-Aj) 


(4.18) 


where Uoo = 1. 

The pressure coefficient for an inviscid and incompressible calculation of M = 0.6 
flow over a 6% thick (r = 0.06) parabolic arc airfoil is shown in Figure 4.3. This 
result was computed solely to check the DP code in an incompressible calculation. 
The linear theory result plotted in Figure 4.3 is obtained for the parabolic arc airfoil 
by a method given in Ashley and Landahl (1965). That method gives the pressure 
coefficient for thin airfoils in incompressible flow. It can be shown that the pressure 
coefficient for a subsonic compressible flow is related to the pressure coefficient for 
incompressible flow by the factor f3 = \J\ — Mjx>: 

G P = \ ( G p ) incompressible (4,19) 

Hence, the DP incompressible flow result in Figure 4.3 is scaled to a compressible 
result using Equation 4.19 just as the linear theory result is scaled from Ashley and 
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Landahl’s incompressible calculation. A compressible stream function solution on the 
same 61 x 41 grid lies right on the scaled incompressible solution for the parabolic 
arc airfoil case. A similar result is shown for a 4% thick (PP = 0.02) sine wave arc 
airfoil in Figure 4.4. The incompressible solution using the dual potential code scales 
to the compressible result of a stream function solution. The computed incompress- 
ible solutions here have been scaled to the compressible result for comparison with 
calculations obtained from available compressible stream function codes. 

4.3. 1.2 Steady viscous flow 


Channel inlet The two-dimensional channel inlet flow has been com- 
puted by many researchers (Wang and Longwell 1964; McDonald et al. 1972; Tenpas 
and Pletcher 1987). This case provides the first test of the vorticity transport equation 
and the vector potential solver. Also, stretched grids must be used in this problem 
for the first time to adequately resolve the inlet features. This problem will be a 
stepping stone to the heat transfer cases presented subsequently. It will be seen that 
the dual potential formulation can solve this flow problem for either a rotational or 
irrotational inlet condition. 

The developing flow in a two-dimensional inlet has been computed for Re = 10 
to 7500. The Reynolds number is based on the hydraulic diameter, 

4 X cross-sectional area 
wetted perimeter 

For a two-dimensional passage the hydraulic diameter is twice the wall spacing, h. 
The physical distances, x and y , for the two-dimensional channel cases are non- 
dimensionalized by the hydraulic diameter. The non-dimensional wall separation 
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distance is then ymax = Note that y for a two-dimensional channel problem will 
then be in the range 0 < y < rj for the results plotted here from wall to wall. The 
geometry and boundary conditions are given in Figure 3.3. A uniform inlet velocity 
profile was prescribed. The inlet conditions can be made to correspond to a rotational 
or irrotational inlet flow. The irrotational inlet may be more realistic for experimental 
models with a rounded entrance (Van Dyke 1970). 

The centerline velocity development is shown in Figure 4.5 for Re = 300. Good 
agreement is obtained for both the irrotational and rotational inlet conditions. Note 
that the x used in the axis label, x/(Re*Dhyd)* 16, in Figure 4.5 is dimensional. This 
was done to form the dimensionless group used by other researchers. The multiplier 
of 16 is used to give the same x axis range as others who have used the channel half 
height instead of the hydraulic diameter in the Reynolds number and dimensionless 
axial length. 

For the results shown in Figure 4.5, all convective terms have been central dif- 
ferenced. Refining the grid for the irrotational inlet case shows that the numerical 
algorithm is clearly better than first order accurate but is not second order accurate 
on a stretched grid. The truncation error is estimated to be approximately O ( Ax)^ 
where Ax represents the x and y grid spacing here. (The solution for the 41 x 41 
grid is approximately one-third of the way vertically between the 21 x 21 grid solution 
and the 81 x 81 grid solution). The need to upwind difference the convective terms 
becomes necessary as the Reynolds number increases in order to maintain stability. 

The skin-friction development is shown in Figure 4.6 for the rotational and ir- 
rotational inlet conditions. Both computational results asymptote to the expected 
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fully- developed Cy Re = 24.0 for the two-dimensional channel. 

The complete Navier-Stokes equations for incompressible flow were used for these 
calculations. It is permissible to drop the streamwise second derivatives for high 
Reynolds number flows. Dropping the u> X x terms causes the flow to develop slightly 
faster (i.e., in a shorter x distance). 

Constant property heat transfer Heat transfer cases have been 
computed for constant properties. The results are in good agreement with published 
computational results (Schade and McEligot 1971; Hwang and Fan 1964). Both 
references above obtained solutions for this problem utilizing boundary-layer approx- 
imations. 

The heat transfer cases can be classified as in Figure 4.7. As shown in the figure, 
there are three basic types of problems for specified wall temperature and/or wall heat 
flux for the parallel plates geometry. Test cases will be computed for wall conditions 
that are constant with x. The most detailed presentation will be for the constant 
wall temperature case. 

Any combination of the wall boundary conditions in Figure 4.7 can be solved by 
the dual potential code. The boundary conditions may be functions of x. 

The notation used in this heat transfer section is the same as used by Shah 
and London (1978) and Kays and Crawford (1980). The x distance in the following 
heat transfer results is referred to the Peclet number as in Shah and London (1978). 
Therefore, 

x* = ^~, Pe = RePr 

A quadratic fit to the computed temperature profile was used to compute the local 



Centerline Velocity Development 
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Figure 4.5: Centerline velocity development for a 2-D channel inlet 








Figure 4.7: Basic wall heat transfer boundary conditions 


Nusselt number, where 


Nu = 


(4.20) 


The Nusselt number was computed from the temperature solution. The notation of 
Kays and Crawford (1980) will be used in the following derivation. Specifically, an 
overdot indicates a time derivative, the subscript o will mean at the wall and a double 


prime ( ft ) will indicate per unit area. The wall heat flux can then be written as 


% - M^wall - ^mean) 

4 ? = -k~ 

°y wall 


(4.21) 

(4.22) 


The wall heat flux is given by either the convection heat transfer equation or the wall 
conduction equation as written above. These two must be equal, so 


M^wall “ Tmean) = -k- 


(4.23) 


The definition of Nusselt number from Equation 4.20 can then easily be written in 
terms of the computed temperatures by rearranging Equation 4.23. The result is 


? hyd = Mgall hyd 
k (^wall ~ ^mean) 


(4.24) 
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where T me an is the local mixed mean fluid temperature. It is also sometimes referred 
to as the bulk fluid temperature or the mixing cup temperature. T m ean is computed 
as the single (mean) temperature for the axial convection that equals the integrated 
axial convected thermal energy rate. That is, 


mcpTmean = {pA c V) cpT m e&n 

= Ja pucpT dAc 


(4.25) 

(4.26) 


where m is the mass flow rate, A c is the cross-sectional area at the x position for 
which Tmean is being computed and V is the mean velocity. Therefore, 

Sa c uT dAc 


Tmean — 


Sa c u <i4c 


(4.27) 


The constant property cases to be presented have fully- developed temperature 
solutions which are self-similar. The self-similar temperature function is 


0 _ r wall Hi lH (4.28) 

^wall — ^mean 

The variable, 0, will be referred to as the temperature parameter in the following 
results. It is defined and used in the dual potential code to indicate when the temper- 
ature solution has become fully developed. This is in addition to an error tolerance 
on the temperature at each point. 

It is important to pack many points close to the channel inlet, so the grid was 
stretched until grid independent solutions were obtained for a 41 x 41 grid. Those 
stretching parameters were then used for further grid refinement by adding points. 
The length of channel was also chosen so as not to interfere with the natural flow 
development. The fully-developed skin friction for these two-dimensional cases is 
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G fRe = 24.0 as computed analytically and by the dual potential code. It is inter- 
esting to note that all of the results for the heat transfer cases were computed for 
the rotational inlet condition. This is the inlet condition used by others reporting 
computational results. Van Dyke (1970) reports that the irrotational inlet condition 
may be more realistic for comparison with experimental models having a rounded 
entrance. The dual potential code can easily handle either a rotational or irrotational 
inlet condition. 


Constant wall temperature Some comparisons of con- 
stant property heat transfer results for constant and equal wall temperatures are 
shown in Figures 4.8 and 4.9. The published results of Schade and McEligot (1971) 
are represented in each figure as a solid line. The present results using the dual 
potential code are indicated by symbols. The grid size for the dual potential results is 
indicated in the figure legend. The Reynolds number based on the hydraulic diameter 
is 150 for the present results in Figure 4.8 and 300 for the present results in Figure 4.9. 
Schade and McEligot use the boundary-layer assumptions and drop their numerical 
results for approximately the first 20 points in x. The dual potential results for the 
161 x 81 grid have only the first 4 points dropped. There is nothing wrong with those 
4 points, they are just well ahead of the data presented in the literature. The dual 
potential incompressible code makes no assumptions other than constant properties. 
This accounts for the Nusselt number discrepancy between the two solutions near the 
inlet. The dual potential code predicts a higher heat transfer rate near the inlet than 
the results of Schade and McEligot. Figures 4.8 and 4.9 show that there is a Peclet 
number dependence for the Nusselt number. For the range of results presented, a 
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Figure 4.8: Local Nusselt number for flow between parallel plates at constant temperature. 
Pe=105. Asymptotic Nu=7.541 











Parallel Plates Constant Temperature 

rotational Inlet, Re-300 












102 


Peclet number of about 210 is required to give good agreement with the boundary- 
layer assumptions. Lower Peclet numbers give a higher heat transfer rate near the 
inlet. Profiles of the temperature parameter, 9, are plotted in Figure 4.11 for several 
streamwise positions. The fully-developed profile is parabolic. 

To match the boundary- layer assumptions used by Schade and McEligot, the 
streamwise second derivatives must be dropped. The affect is similar to assuming a 
high Reynolds number flow. Schade and McEligot report that Pe > 100 is sufficient 
for the boundary-layer assumptions to be valid. Also, Schade and McEligot neglect 
viscous dissipation. For incompressible flow, the only streamwise second derivatives 
in the dual potential method are u>xxi Axx and Txx- When these derivatives and the 
viscous dissipation are dropped, the results approach those of Schade and McEligot, 
as can be seen in Figure 4.10. It is clear that the boundary-layer assumptions cannot 
be used to get accurate results near the inlet at lower Reynolds numbers or Peclet 
numbers. The boundary-layer assumptions neglect the important elliptic effects at 
the inlet and underpredict the heat transfer. 

The computation rate for the dual potential code is 80-100 MFLOPS overall 
using a single processor on a Cray X-MP. The computation rate can be stated in a 

more useful manner as ~ global iteration^ grid point on a Single P rocessor of the 
Cray X-MP. The number of global iterations is on the order of the number of grid 

points. A conservative estimate of Cray X-MP cpu time is: 

• 2 

cpu time (s) = 5.5 x 10~° x (number of grid points) 

The two-dimensional dual potential code is a compressible code. The compressible 
terms are computed for this calculation even though they are zero. If the code was 
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designed purely for incompressible solutions the code would be much faster — at least 
twice as fast based on a rough operation count. The error tolerance used for these 
calculations was e < 10 — 

The solution to the constant property case is independent of wall temperature. 
This was checked using T wa jj = 2.0,1.1,1.01,1.001,1.0001 with the same constant 
property results obtained for each case. Shah and London (1978) also present data 
for the constant wall temperature case. For the present results, the grid was stretched 
using <r = 1.05, a = 0.5 and (3 = 1.2. 

Constant heat flux The constant temperature case has 
revealed that Pe = 210 gives good agreement with the boundary-layer approximations 
used by Schade and McEligot. Just one case will be presented here in Figure 4.12. The 
results compare well to solutions obtained using the boundary-layer approximations, 
but are more accurate near the inlet. A second order polynomial curve fit to the 
temperature distribution was used to obtain the wall temperature. The temperature 
parameter profiles for this case are presented in Figure 4.13. 

Mixed wall boundary conditions Solutions were com- 
puted for two cases with one wall insulated and the other wall at either constant 
heat flux or constant temperature. Shah and London (1978) refer to these boundary 
conditions as the fundamental boundary conditions of the second kind and third kind 
respectively. The asymptotic Nusselt number for the case of one wall insulated and 
the other at constant heat flux is computed analytically to be Nu = 5.385. This case 
will be designated by H2 in the results. The asymptotic Nusselt number for the case 
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Figure 4.10: Local Nusselt number for flow between parallel plates at constant tern 

perature. Dual potential code with boundary-layer assumptions 
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Figure 4.11: Profiles of the temperature parameter at various 

wall temperature case 
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Figure 4.13: Profiles of the temperature parameter at various s* for the constant 

wall heat flux case 
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Figure 4.14: Local Nusselt number for flow between parallel plates with mixed wall conditions. Pe=210 
Asymptotic Nu=5.385 for one wall insulated and the other at constant heat flux. 
Asymptotic Nu=4.861 for one wall insulated and the other at constant temperature 




Temperature Parameter Profiles 

rotational Inlet, Re^OO, H2 wait conditions 



Figure 4.15: Profiles of the temperature parameter at various z* for the case of one 

wall insulated and the other at constant heat flux 
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of one wall insulated and the other at constant temperature is computed analytically 
to be Nu = 4.861. This case will be designated by H3 in the results. The computed 
Nusselt number development for these two cases and the analytical asymptotic Nus- 
selt number are shown in Figure 4.14. Temperature parameter profiles for these two 
cases are plotted in Figures 4.16 and 4.15. 

Boundary layer This case is considered to provide yet another check 
on the accuracy of the dual potential scheme since the laminar boundary-layer solu- 
tion is well known. Despite the simplicity of this flow, it is not always easy to obtain 
good comparisons with numerical solutions to the full Navier-Stokes equations. 

The boundary conditions for the laminar incompressible flow over a flat plate are 
shown in Figure 4.17. Since this is an incompressible constant property flow, B,T 
and p need not be computed. 

Concerning the vorticity boundary conditions, an irrotational or rotational inlet 
flow may be specified. The best comparison with Blasius results was obtained for 
a zero inlet plane vorticity (irrotational condition). This gives a non- zero v velocity 
at the inlet. Both u and v are fixed to be zero at the stagnation point. The wall 
vorticity is computed as in all viscous wall cases. The freestream is assumed to be 
irrotational. The exit condition is based on the usual fully-developed assumption that 
second order streamwise derivatives are zero. 

The vector potential inlet and viscous wall conditions are the usual A = 0. In this 
case the vector potential is used to account for the throughflow velocity (i.e., leakage 
velocity) at the top boundary. It was stated in the boundary condition section that 
only the scalar potential would be used to give throughflow velocities. This is the only 



Figure 4.17: Boundary conditions 
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= 0 

A X = 4>y — V 
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U = 1 
Vy=0 


exit: 
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Uxx = 0 
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or laminar incompressible flow over a flat plate 
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exception. The v velocity is given at the top boundary by assigning the same velocity 
at the boundary as at the point just inside. This implies that ^ = 0. Hence, Ax 
along the top boundary is computed from the assumed v leakage velocity using the 
velocity decomposition for v, v = <f>y — Ax- More elaborate calculations can be used 
to compute the v leakage velocity from the displacement thickness. For example, the 
leakage can be accurately written as 

dS* v 

dx Uoo 

Thus the leakage can be controlled based on a known or computed displacement thick- 
ness, S*. These ideas worked no better than simply assigning v(i,nj) = v(i,nj - 1) 
at the freestream boundary. At the exit u* = 0 was acceptable, but uxx = 0 gave 
a slightly better boundary-layer profile as compared to the Blasius (1908) solution. 
Note that the exit condition on the vector potential, Axx = 0, is equivalent to an 
extrapolation technique. The vector potential at a point just outside the exit plane 
boundary is computed by Axx = 0 for use in the Poisson equation for A. The Poisson 
equation is solved at the exit plane. This assumes that the flow is fully developed and 
removes the need to specify the velocity at the exit. The scalar potential boundary 
conditions are the same as for the incompressible channel case. Therefore, the solu- 
tion of the scalar potential Poisson equation is known a priori and <j> ■= x is simply 
assigned and not re-computed. The same boundary conditions as the channel case 
were used to avoid a solution for the scalar potential in this case. The only diffi- 
culty was a need to change the top boundary condition on the vector potential to 
account for a leakage flow. This leakage flow is small and could be neglected with the 
boundary far enough away. If it is neglected and the top boundary is too close, the 
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boundary-layer profile tends to overshoot the freestream value by several percentage 
points. Allowing leakage at the top wall permits the solution to be obtained on a 
smaller grid. 

The laminar boundary-layer computation agrees well with the Blasius solution. 
A boundary-layer profile comparison is shown in Figure 4.18. An 81x81 grid was 
used with stretching in the y direction. This puts about 22 points in the boundary 
layer at Rex = 312,500. The solution is attained in 54.5 seconds at 85 MFLOPS 
overall on a single processor of the Cray X-MP. The data marked with triangles 
were computed by limiting the solution domain of the vorticity transport equation to 
31 transverse points. This example was specifically used to demonstrate the ability 
to reduce the computational effort by dropping the calculation of vorticity where 
it is negligible and still obtaining a full field solution for the problem. Even faster 
solutions are possible using this method. Normally the Poisson equations dominate 
the solution time. That does occur in this problem also, but iterations on the vector 
potential calculation were limited with no deleterious effect on the final solution. With 
the vector potential limited, the vorticity transport equation becomes the dominant 
time consuming computation for this problem. It is this condition that was used to 
demonstrate the ability to reduce computation time for the vorticity. This speed-up 
affect is shown in Table 4.1. 

For problems that one knows have vorticity approaching zero at a boundary, 
the point-wise convergence check used in the code can be helpful in deciding where 
and how much to limit the solution domain for vorticity. For example, if vorticity 
goes to zero at the freestream boundary as in this problem, one will note that for the 
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Table 4.1: Speed-up affect of limiting the vorticity solution domain 



Transverse pts for vorticity solution 

CPU time (s) 

MFLOP 

81 

54.5 

85 

41 

40.2 

77 

31 

36.9 

73 


convergence checking used here (which is point by point and normalized on the largest 
value in a plane), the biggest error will be in the region where vorticity is important. 
If the point of maximum error is far away from the known uj = 0 boundary, then 
limiting of the solution domain is indicated. 

The MFLOP rate in Table 4.1 drops as the solution domain for vorticity is limited 
because the vorticity transport subroutine had a high MFLOP rate. The solution was 
speeded up by limiting the vorticity solution domain and thereby reducing the CPU 
time spent in computing the vorticity in regions where it is nearly zero anyway. The 
vorticity transport subroutine is a 107+ MFLOP code. Reducing the computation 
time on vorticity reduced its high MFLOP rate in the weighted average to compute 
the overall MFLOP rate for the code. 

Solutions were obtained for Re® up to 2,000,000 using both rotational and irrota- 
tional inlet conditions. The rotational inlet condition gave a higher skin friction and 
a slight velocity overshoot. The computed skin friction along the flat plate for the 
two inlet conditions is shown in Figure 4.19. The rotational inlet condition produces 
a higher C j than the irrotational inlet by the conservation of vorticity. Since the 
inlet vorticity is non-zero in the rotational case, there is an additional amount of vor- 
ticity in the rotational as opposed to the irrotational inlet case. The boundary-layer 
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calculation for the rotational inlet (u> ^ 0) has vorticity in excess of that assumed in 
the Blasius solution. The extra vorticity is the inlet amount. The Blasius solution 
is obviously for an irrotational inlet. The results presented here were obtained using 
first order upwinding of the streamwise convective term. Weighting toward central 
differencing works but requires more CPU time to converge. The equations used here 
do not use the boundary-layer assumptions and so can accurately model the near 
leading edge flow. The skin friction is computed as: 


du 

/ ^ wall ^7 

C f yJ Re x = — U 


wall / Poo^oo x 
Moo 


(4.29) 


ypo o Voo ' 

The velocity derivative is computed using a second order polynomial fit to the u ve- 
locity component. The computed asymptotic value for the skin friction is C jyj Rex = 


0.664 for the irrotational inlet case and C j V^a: = 0.704 for the rotational inlet. 

The x and y distances in Figures 4.19-4.21 are non-dimensionalized by the length 
of plate required to give Re^ = 312, 500. Velocity vectors for this boundary-layer flow 
are plotted in Figures 4.20 and 4.21. The most easily observed difference in these two 
figures is the flow at the inlet plane. The irrotational inlet flow condition produces a 
v velocity component at the inlet plane. 





Skin friction along a flat plate, ReL=312 500 
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Figure 4.19: Skin friction along a flat plat 
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4.3.2 Compressible flow 

4. 3. 2.1 Steady irrotational inviscid flow 

Bump cases The flow assumptions for this problem are the same as 
in the previous bump test case: 

1. inviscid 

2. irrotational 

3. isentropic 

4. no body forces 

The flow can be incompressible or compressible and steady or unsteady. The bound- 
ary conditions for this case are shown in Figure 4.22. The same boundary conditions 
apply for the unsteady problem which is discussed in the next section. The assump- 
tions for this test case are the same as in the bump cases previously discussed on 
pages 83-92 except that the compressible solution is computed directly rather than 
scaled from an incompressible result. 

There are many possible well-tested solution methods for this irrotational prob- 
lem. The dual potential method is easily adapted to solve the flow using any of the 
following combinations of governing equations. Either potential may be used. The 
scalar potential may be used as the traditional velocity potential for a full potential 
solution. The vector potential may assume the role of the stream function for a com- 
pressible stream function solution. In addition, an Euler mode of solution is available 
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5 = vy initially, then by continuity 

y,j index 

freestream: 

p from Bernoulli 
T from constant total enthalpy 

u = 4> x 



PslH 

II 




bump 


inlet: 
5 = 0 

<j> = 0 
P = 1 
T = 1 
u = 1 
v = 0 



exit: 

5 = 0 
<^1 = 1 
P = 1 
T = 1 
u = 1 

V = 4>y 


freestream: 

B from continuity index 

4>y - 0 

p from Bernoulli 
T from constant total enthalpy 

u = <f>x 
v = 0 


Figure 4.22: Boundary conditions for laminar compressible irrotational flow over a 

bump 
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using the dual potential method with viscous terms and thermal conductivity set to 
zero. 

The traditional potential methods simply require the velocity potential or stream 
function with a version of the Bernoulli equation. These two methods will be pre- 
sented in a form only applicable to steady flows. The full potential equation written 
in conservation law form is 

(^*) + ^ ( Ph ) = 0 ( 4 - 30 ) 

This may be solved by the scalar potential solver in the dual potential code with a 
redefinition of the scalar potential derivatives to be the parenthesized quantities in 
the equation above. This requires a modification of the boundary conditions also. 
Since the flows to be computed here are subsonic, it is permissible to write the full 
potential equation in non-conservative form: 

P 

This is easily solved by the scalar potential solver with no modification to the bound- 
ary conditions. The steady compressible Bernoulli equation is then solved for density, 

+ ( 4 . 32) 

Using the stream function, only the velocity decomposition and Bernoulli equa- 
tion are required. For the stream function approach, solve the following: 
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T = T t - 


2 , 2 
u* + v* 

2cd 


The stream function, ip, above is defined to satisfy continuity automatically by: 


u = 


4>y 


_ ipx 
P 

The above would require some code modification to define the vector potential, 
A, as the stream function, ip. An alternate approach is to define 

A _ ** 

A-X — 

P 

Ay = — 

^ P 

then the vector potential solver remains unchanged for a stream function solution and 
only the boundary conditions must be modified by p. This also permits easy use of 
the non-linear form of Bernoulli’s equation which becomes 

7 - 1 


= {i + (1-^-4)}''- 


(4.33) 


The energy equation is not specifically needed if the Bernoulli equation is solved since 
it is a statement of conservation of momentum and energy. 

The dual potential formulation may be easily converted to an Euler solver by 
dropping the viscous terms and thermal conductivity. All that is required is to set 
p — k = 0. For the bump problems, the flow is irrotational so that the vorticity 
transport equation is also dropped. The DP code in the Euler mode solves the 
governing equations in the following order: 


1. V 2 <P = B 
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2. V = V 4> 

3. the energy equation for T 


4. an equation for B 

5. an equation for p 

6. the equation of state, p = pRT 

The assumptions for this problem allow the energy equation to simplify to a 
statement of constant total enthalpy. There are several equations to choose from to 
compute density and dilatation. Using the DP code in Euler mode, three possible 
equation sets are given below. These are listed in order of increasing complexity, or, 
in other words, from fastest to slowest. 

p by Bernoulli 
/ B from continuity 


p from continuity 
B by its transport equation 


p by Bernoulli 
B by its transport equation 


U = V <j> 
V 2 4> = B 


All these equation sets give the same results. Obtaining p from the Bernoulli 
equation and B from continuity is as fast as the compressible stream function ap- 
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proach. Using the continuity equation to obtain p and the dilatation transport equa- 
tion for B takes about 10 times more CPU time than the fastest methods. Only the 
last equation set will be applicable to an unsteady flow, provided the energy equation 
is also time accurate. 

For the compressible flow over a bump, the dual potential code uses the following 
initial conditions: 

1. u = 1, v = 0 except at the bump where v = & (the slope of the boundary) 

2. B = 0 everywhere except at the bump, where B = vy (Euler version only) 

3. p — 1 everywhere 

4. <j> = 0 everywhere (Euler and full potential versions only) 

5. A = 0 everywhere (Stream function version only) 

6. T = 1 everywhere 

7. p = k = 0 everywhere for all time. 

There were initial difficulties solving the bump test cases due to abrupt function 
changes at the boundaries. It was necessary to smooth the variables ( T,p,B ) into 
the boundaries either by extrapolation or by solving for the variable at the boundary 
using available information and a governing equation — usually a simplified equation 
such as the Bernoulli equation. The far- and near-field boundaries then have values 
that are near the expected steady condition but differ slightly to provide smooth 
derivatives into the boundary. For example, the density is not exactly 1.0 at the 
far-field boundary when it is computed by the Bernoulli equation at that boundary, 
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but the derivatives of density which occur in the dilatation transport equation and 
elsewhere are much smoother using this approach. The y=constant boundaries were 
the most critical. 

The pressure coefficient can be computed from the definition 

n P-Poo 

P = K 1/2 

%Poo too 

or, if the scalar potential is used, 

C p = - 2 <j>i - 2 <t> x - <f>y 

where + v^. 

The flow over a 4% thick sine wave arc airfoil and a 6% thick parabolic arc airfoil 
was computed for a freestream Mach number of Moo = 0.6. The dual potential code 
was used in a compressible stream function mode and in an Euler mode. Results 
for the sine wave arc are shown in Figures 4.23-4.28 and results for the parabolic 
arc in Figures 4.29-4.34. The DP code is solved on a uniform 61 x 41 grid. The 
computational domain is 3 chords in the streamwise direction and extends 2 chords 
away from the airfoil. The airfoil is centered in x. The sine wave arc was the first 
successful compressible calculation. There are no stagnation points to contend with 
on the sine wave arc. Note the fair agreement between the stream function and Euler 
solution for this case shown in Figure 4.23. A grid refinement study shows that the 
61 x 41 grid gives an Euler solution independent of further refinement. The solutions 
have not been obtained on a larger domain to examine whether the solution is grid 
independent in that respect. The discrepancy between the stream function result and 
the Euler solution is caused in part by the inconsistency of the small perturbation 
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boundary condition with the stream function and Euler solvers. The intent here is to 
show the capability of the solver, not ultimate solutions. The Euler solution requires 
about 3-10 times more CPU time than the potential solution using the same DP 
cot l e — depending on which equations are used to solve for the density and dilatation. 
In the Euler code the convective terms are weighted toward central differencing by 
0.9. The Euler solution is not exactly symmetric as can be seen in the contour plots 
of dilatation (Figures 4.24 and 4.30). This has been observed by others for inviscid 

solutions. 

4. 3. 2. 2 Unsteady irrotational inviscid flow 

Bump cases The flow assumptions for this problem are the same 
as in the previous bump test case, except that the flow will be unsteady. The dual 
potential code will be used in an Euler mode. Constant total enthalpy has been 
assumed for the energy equation. Formally this requires a steady flow, but it will be 
applied here at each time step for the unsteady problem. The particular unsteady 
problems to be computed here have primarily low frequency disturbances (Ballhaus 
and Steger 1975; Beam and Ballhaus 1975). The validity of using the energy equation 
in a quasi-steady fashion was evaluated by computing the unsteady flow with the 
complete energy equation. The pressure coefficient results were within 4%. This 
problem was studied mainly to test the dilatation transport equations and the density 
determination from the continuity equation. The constant total enthalpy (or constant 
total temperature since cp is constant) is written 


(4.34) 



.6 FLOW OVER R THIN RIRFOIL 


129 



sro 


SO'O 



o 


OZ'O 


oro 


— r~ 

oo - o 


so'o- oro- 


Figure 4.23: Pressure coefficient for Moo = 0.6 flow over a sine wave arc airfoil with 

maximum thickness = 4% of chord 
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Figure 4.27: V velocity profiles from inlet to exit for the sine wave arc airfoil 
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Figure 4.28: Comparison of the dilatation profiles for the stream function and Euler mode 
of computation by the DP code. Sine wave arc airfoil 
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Figure 4.29: Pressure coefficient for Moo = 0.6 flow over a parabolic arc airfoil with 

maximum thickness = 6% of chord 
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Figure 4.31: Density profiles from inlet to exit for the parabolic arc airfoil 





Dilatation at the top wall 

parabolic arc, height - 0.03 
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Figure 4.32: Dilatation profiles from inlet to exit for the parabolic arc airfoil 




139 


o 



I 1 1 1 1 1 1-7 

sro oro so - o oo - o so'o- oro- sro- 




Figure 4.33: V velocity profiles from inlet to exit for the parabolic arc airfoil 
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Figure 4.34: Comparison of the dilatation profiles for the stream function and Euler mode 
of computation by the DP code. Parabolic arc airfoil 
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The unsteady case computed here is that of an airfoil thickening from zero thick- 
ness to 10% of chord. The airfoil grows in 20 chord lengths of time measured with 
respect to the freestream velocity. The airfoil growth is given in Figure 4.35. For this 
problem r m ax = 0.1 and t 0 = 20. This test case may be used as a model for the flow 
about the advancing rotor of a helicopter. The dual potential results are compared 
to an Euler solution and transonic small disturbance equation solution of Beam and 
Ballhaus (1975) in Figure 4.36. The pressure coefficient as a function of time is plot- 
ted at the position x/c — 0.525 measured from the leading edge of the airfoil. The 
dual potential code in Euler mode computes the unsteady solution for the thickening 
parabolic arc airfoil in best agreement with the small disturbance computation. All 
results in the figure are for a domain that is three chord lengths long and extending 
two chord lengths out from the bump. The grids are nearly alike. The DP code had 
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61 streamwise by 41 transverse points all uniformly spaced. With a 121 x 81 grid on 
this domain the DP solution is only slightly lower (no more than 1.3%) with a more 
pronounced overshoot. Beam and Ballhaus used 25 points on the airfoil. No claim 
is made that this solution is grid independent. In fact, further calculations using the 
DP method on refined grids and on a domain that is five chords in the streamwise 
and four chords to the far field have given pressure coefficient results that asymptote 
about 8.5% below the results in Figure 4.36. Though the DP Euler solution is short 
of the overshoot of Beam and Ballhaus, it is close to the transonic small disturbance 
(TSD) equation result for this case. 

The solution for this unsteady problem converges slowly at first as the scalar po- 
tential is computed from an initial field of <f> = 0. Later, during periods of slow change, 
the solution is obtained rapidly. When the bump grows the fastest the solution speed 
is slowest. 

Another result is shown in Figure 4.37. For this test the parabolic arc airfoil 
grows to 10% thickness in t = 15 chord lengths. The dual potential solution is 
compared to a solution of the linearized transonic small disturbance equation. The 
pressure coefficient is plotted for x/c = 0.525 as measured from the leading edge 
of the airfoil. Early unsteady calculations made it to steady state by reducing the 
tolerance on the dilatation, B. The boundary conditions for the unsteady problem 
are actually fixed for all time except at the bump. Thus, unsteady terms are not 
required in most boundary conditions at the far- and near-field. A tight tolerance on 
the scalar potential is important here. 
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Figure 4.36: Variation of the pressure coefficient with time at x/c = 0.525 for 

thickening parabolic arc airfoil. Moo = 0.785, maximum thickness 
0.1 in t = 20 chord lengths 
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4. 3. 2. 3 Steady viscous flow 


Variable property channel flows Only a few results will be shown 
to demonstrate the capability of the dual potential code. Clearly, an infinite num- 
ber of different wall boundary conditions exist. The full Navier-Stokes solution for 
variable property flow in a two-dimensional channel could be computed for any wall 
conditions of interest. The computation of variable property channel flows has not 
been wholly successful. High wall temperatures and high wall heat fluxes lead to a 
computed flow separation near the inlet and subsequent reattachment and flow devel- 
opment. The addition of energy to the flow gives a lower “fully-developed” Nusselt 
number than the constant property cases. Also, the Nusselt number does not asymp- 
tote as it does in the constant property case. Rather, the Nusselt number continues 
to drop for wall heating. 

Results for Re = 40 and Re = 150 will be presented for the constant wall 
T 

temperature case with m wa ^ =1.1. The higher Reynolds number is more difficult 

1 inlet 

to compute. Contours of temperature, u velocity and dilatation are shown in Fig- 
ures 4.39, 4.40 and 4.41 for the Re = 40 case. Similar plots for Re = 150 are shown 
in Figures 4.42, 4.43 and 4.44. The Nusselt number development for these cases is 
plotted in Figure 4.45 along with results from Schade and McEligot (1971) and a 
solution from the code of Nelson and Pletcher (1974). Both of these references used 
the boundary-layer equations. 

T 

For a wall temperature of T wa ^ = 2.0 the flow separates near the inlet and re- 

1 inlet 

develops. The separation is evident in the contour plot of u velocity in Figure 4.46. 
This case does not converge. A high wall heat flux will also not yield a converged 
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solution. A non-dimensional heat flux, Q , is defined as 

^hyd 


Q = 


k 0 T 0 


(4.35) 


A near inlet separation develops as shown in Figure 4.47 for Q = 25.0. Lower wall 
heat fluxes do yield converged solutions. A non-dimensional wall heat flux of Q = 0.5 
was computed for Re = 150 and M = 0.1. Contours of the temperature, u velocity and 
dilatation are shown in Figures 4.48, 4.49 and 4.50 respectively. The Nusselt number 
was again below constant property predictions and does not asymptote but continues 
to fall. This behavior is expected for solutions which include viscous dissipation (Shah 
and London 1978). 

The difficulties in solving this problem have been isolated to the density and 
dilatation determination. Let us consider the density first. Density is computed from 
the continuity equation. Since that equation has first order time, x and y derivatives, 
only one condition on time, * and y are needed. In Figure 3.3 for the boundary 
conditions of this problem, however, note that density is constrained by an inlet and 
two wall conditions. An initial condition is also assumed. This is one condition too 
many. The specification of the wall density at both walls seems to be the problem. 
Though many methods of determining the wall density are possible, none seemed to 
be completely satisfactory. The wall density has been determined alternatively from 


1. the x momentum equation solved along the wall, 


2. continuity solved at the line of points just above the wall followed by a statement 
that py — 0, or 

3. density equals a constant, p = y r. 
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Of these, number two gives a constant pressure on a transverse plane, but a realistic 
looking negative streamwise pressure gradient has not been obtained. Mass is con- 
served for all of these approaches. All approaches give the same Nusselt number for 
the conditions computed. So, at least for these low Mach number, and “low” Reynolds 
number flows, the density doesn’t seem to interfere much with the temperature field 
determination. 

Next, consider the dilatation variable, B. The dilatation transport equation 
evidenced early on that it would be difficult to solve. Many approaches were tried: 

1. Divergence theorem constraint on the B field 

2. Under-relaxation 

3. Weighting the B field using the transport equation solution and the definition, 
B = Ux + vy from the velocities. 

4. Smaller time steps. 

A useful check is to make sure that the solution for the dilatation satisfies the 
divergence theorem (also called Green’s theorem or Green’s identity). Only two- 
dimensional compressible problems in Cartesian grids were studied so the following 
derivation will be useful for all the compressible test cases. 

From the Poisson equation for the scalar potential, Equation 2.13: 

V 2 <t> = B (4.36) 


Thus; 


ll D v2 * dA = ll D BiA 


(4.37) 
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Applying the divergence theorem yields: 

<4 - 38) 

IJd b dA = JJp (u x + ft/) dA = (un x + vny ) ds (4.39) 

where D is the domain, dD is the boundary of the domain and n is the outward 
unit normal. The implication of the above is that the scalar potential should be used 
to handle the throughflow velocities at domain boundaries. This is exactly how the 
scalar potential boundaries are imposed (except for the boundary-layer calculations 
where the vector potential is used for the throughflow velocity). A representative 
two-dimensional geometry is shown in Figure 4.38. The outward normal derivative is 
indicated on the figure. Proceeding counterclockwise around the domain and applying 
the divergence theorem, the area integral of B must satisfy: 

fj D BiA = -r-**r«* 

4- f —vdx+ [ udx (4.40) 

Jxi Jy\ 

Writing the above integrals in the positive coordinate directions gives: 



The area integral of dilatation and the boundary integrals of the throughflow velocities 
are easily computed to provide a check on the dilatation field. 

This constraint was observed to be satisfied automatically to acceptable tolerance 
for all the compressible problems except the viscous variable property channel flows. 
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Figure 4.38: 2-D solution domain for divergence theorem application 
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In that case, converged solutions can best be obtained by enforcing the divergence 
theorem constraint on the computed B field. 

In spite of these problems, solutions were obtained without applying any sim- 
plifying assumptions (i.e., no terms were dropped in the governing equations). The 
results here were obtained by method 2 for the density boundary condition and the 
divergence theorem constraint applied to the solution for B . 

Some trends that were observed will be reported next. For variable property 
flow, the Nusselt number does not asymptote but continues to fall for wall heating. 
This is in agreement with the findings of others. Shah and London (1978) indicate 
that the limiting Nusselt number is zero for variable property flow with wall heating. 
This is reason for confidence in the trends computed by the dual potential method. 
Dilatation has no choice but to conform to the numerical satisfaction of Green’s 
theorem. This constraint must be used throughout the calculation. The correction 
becomes smaller and smaller as the steady state is approached, but without the 
correction the convergence is slowed and seemingly halted. At this point the only 
remaining problem appears to be getting a believable negative streamwise pressure 
gradient. The finest grid solution presented here is for 41 x 41 points. For this coarse 
grid the skin friction asymptotes to Cf Re = 23.6. The fully- developed skin friction 
for an incompressible channel case would be C yRe = 24.0. 

The code runs at 83 MFLOP for full Navier-Stokes calculations. The low wall 

heating cases converge at - glob al iteration x grid point ‘ 
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Constant wall temperature Converged solutions were ob- 

m 

tained for = 1.1. Results for Re = 40 and Re = 150 will be presented. The 

r inlet 

flow Mach number is M = 0.1. Contours of temperature, u velocity and dilatation are 
shown in Figures 4.39, 4.40 and 4.41 for the Re = 40 case. Similar plots for Re = 150 
are shown in Figures 4.42, 4.43 and 4.44. Recall that the physical distances, x and 
y, for the two-dimensional channel cases are non-dimensionalized by the hydraulic 
diameter. The non-dimensional wall separation distance is then t/max = so that 
y is in the range 0 < y < ^ for the results plotted here from wall to wall. The u 
velocity contours in Figures 4.40 and 4.43 are similar to results for an incompress- 
ible channel case. Notice that the core velocity attains u = 1.5 which is the exact 

value for the incompressible case. The temperature contours in Figures 4.39 and 

T 

4.42 are realistic considering that the wall temperature is a constant m ffa ^ = 1.1. 

1 inlet 

The dilatation contours indicate that most of the compressible effects for this flow 
are concentrated near the inlet. The Nusselt number development for these cases is 
plotted in Figure 4.45 along with results from Schade and McEligot (1971) and a 
solution from the code of Nelson and Pletcher (1974). Schade and McEligot used the 
boundary-layer assumptions and neglect viscous dissipation in the energy equation. 
They compute an increase in Nusselt number for heating when properties are variable. 
The computations of Pletcher are also based on the boundary-layer equations but in- 
clude the dissipation terms in the energy equation. His results confirm the behavior 
predicted by the dual potential solution of the Navier-Stokes equations. Shah and 
London (1978) also confirm that the Nusselt number does not asymptote for variable 
property conditions with heating. 
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Parallel Plates Constant Temperature 

variable properties, rotational Inlet, M-0.1 
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Figure 4.45: Nusselt number comparison for constant wall temperature and variable 
properties 
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Figure 4.46 is from the results of the case with = 2.0. This case does not 

1 inlet 

converge. These results are shown mainly to show the separation that is computed 

by the DP code. For convenience, the constant wall temperature is given as T wa jj in 

T 

the figures. It is to be understood that it is actually the ratio of 

1 inlet 


Constant wall heat flux Figure 4.47 is for the case with 
Q = 25.0. This case also does not converge. These results are shown mainly to show 
the separation that is computed by the DP code. The high heat flux and high wall 
temperature cases exhibit this same feature. 

A converged solution is obtained for lower heat fluxes. Results are plotted in 
Figures 4.48-4.50 for Q = 0.5, Re = 150 and M = 0.1. The non-dimensional wall 
separation distance is ymax = \ since y is non-dimensionalized on the hydraulic 
diameter which is simply twice the wall spacing. 


Compressible boundary layer 

Subsonic freestream The compressible subsonic flow over 
a flat plate was computed to test the code on an external, viscous, compressible flow 
case. The conditions at the plate are adiabatic. The normal pressure gradient at the 
plate, py, is assumed to be zero. This is consistent with a boundary-layer assumption. 
The boundary conditions for this case are shown in Figure 4.51. 

The full dual potential method and a variation were used to compute the solution. 
The variation on the dual potential method was used as a self check on this problem. 
It uses dilatation computed from the continuity equation and density computed from 
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freestream: 


inlet: 

u> = 0 
B = 0 
T = 1 
P = 1 
>1 = 0 
<£ = 0 
u = 1 
v ^ 0 


u> = 0 

B _ ~(tt px+vpy) 

T = 1 
P = 1 

■ 4 * = <f>y — V 

4>y = 0 

u = 1 

Vy = 0 



exit: 
w®® = 0 

£®x = 0 

r®® = 0 

px upwinded 
Axx = 0 
<f> x = u 
Uxx = 0 
Vxx = 0 


B = 0 
dT _ n 

W“° 

flat plate: p — fixed by py = 0 and T distribution 

>1 = 0 

4>y = 0 
u = 0 
v — 0 


x, i index 


Figure 4.51: Boundary conditions for laminar compressible flow over a flat plate 
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Figure 4.52: Boundary layer profile at Re x - 100, 500 and Moo 
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Figure 4.53: Temperature distribution for adiabatic flow over a flat plate at M 
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the ideal gas law as follows. From continuity, and assuming steady state conditions, 

n ~ ( U PX + Vpy} 

B = 

p 

It is possible to assume that the pressure is constant and uniform for this case, so 
that density can be computed from the ideal gas law by 


P = 


P 

RT 


where = 1 

T R 


This variant is possible because of the simplifying assumptions of this problem. This 
gives a much faster solution than the transport equations. When the transport equa- 
tions are solved for dilatation and density, the boundary conditions shown in Fig- 
ure 4.51 are used. The wall density can actually be obtained by other methods as 
mentioned in the variable property channel flow discussion on page 146. A compar- 
ison of the two solution methods and results from a boundary-layer finite-difference 
scheme (Christoph and Pletcher 1983) are shown in Figure 4.52 for a M = 0.5 flow. 
Pletcher’s data were computed using the above mentioned finite- difference boundary- 
layer scheme with approximately 50 points in the boundary layer. For reference, the 
Blasius profile and a solution at M = 1.0 from Schlichting (1979) are included in Fig- 
ure 4.52. It can be seen that the compressible boundary layer thickens as the Mach 
number increases. 

The temperature profile for this case is shown in Figure 4.53. The results from 
the boundary-layer finite- difference scheme are not plotted but coincide with the 
Pr = 0.7 theoretical curve. More grid points in the dual potential solution could 
be required for better agreement. However, the results are within 1% of theory. In 
the present calculations, the wall temperature is computed by using the zero wall 
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derivative condition for temperature as a boundary condition on the energy equation. 
Pletcher reported setting the adiabatic wall temperature in his calculations. The 
wall temperature is 1.042 by theory and 1.039 as computed using the dual potential 
method. 

This was the first compressible viscous case to be computed successfully. There 
was initially a problem computing this flow using the full dual potential equation 
set. The dilatation transport equation could only be converged by dropping the 
terms which originated from the pressure terms in the primitive variable momentum 
equations. (Dropping those terms is equivalent to assuming that the pressure is 
everywhere constant.) This led to the determination that those terms (second order 
derivatives in the B equation source term) must be differenced conservatively. This 
is necessary because the pressure gradient is zero, yet numerically it is non-zero for 
p = pRT substituted into the pressure derivative terms. 

These cases were computed by the full dual potential method without using 
boundary-layer assumptions other than py = 0 at the plate. The full energy equation 
was used. Streamwise convective terms were upwinded. The time step limitation is 
controlled by the dilatation variable. The skin-friction development is nearly indis- 
tinguishable from the incompressible boundary layer solution for an irrotational inlet. 
The skin friction asymptotes to a value slightly lower than for incompressible flow 
with the difference in the fourth place behind the decimal point (c/. Van Driest 1952). 
The asymptotic value computed here was Cf\/ R^ = 0.6646 at Re x = 100,500. At 
the same Reynolds number the incompressible solution computed C f %/Re = 0.6647. 

With 81 transverse points it is observed that the vorticity magnitude is less than 
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10 -5 from j = 35-81 and the dilatation magnitude is less than 10 — ^ from j = 54-81. 
The solvers for vorticity and dilatation can be shut off beyond these points with no 
adverse affect on the solution and a speed up of approximately 14%. 

Supersonic freestream A supersonic flow over a flat plate 
was also computed. The conditions at the plate are adiabatic. The result in Fig- 
ure 4.54 was obtained by computing the dilatation from the continuity equation. A 
zero streamwise pressure gradient was assumed. Then, by the boundary-layer assump- 
tion, a zero normal pressure gradient is also assumed so that pressure is constant for 
this problem. The dilatation may then be computed from the continuity equation 

B = (4.42) 

P 

and density can be computed from the ideal gas law and the constant pressure as- 
sumption 

P 1 

p = -i— = — for a unit pressure field. (4.43) 

RT T 

It was not possible to compute this flow using the transport equation for B. 

The abscissa in Figure 4.54 is where ^ is the momentum thickness defined 

2 

Experimental results for this case were obtained by O’Donnell (1954). The theoretical 
calculation is by Chapman and Rubesin (1949). The achievement of good results here 
may be due in part to the use of the fixed pressure assumption. 
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Figure 4.54: Velocity profile in an adiabatic laminar boundary layer in supersonic flow 
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Figure 4.55: 3-D duct geometry 

4.4 Three-dimensional Cases 

Only a first step to three-dimensional calculations using the dual potential method 
has been taken here. For incompressible flow, solutions were computed for the devel- 
oping flow in rectangular ducts of constant cross section. This may be considered as 
a stepping stone to more complex geometries. Compressible flow solutions in three- 
dimensions have not been obtained. Time did not permit further development of the 
dual potential method, but the capability of this method has been demonstrated for 
three-dimensional incompressible flow calculations in a simple geometry. 

4.4.1 Incompressible flow 

4.4. 1.1 Incompressible channel inlet A three-dimensional, laminar, in- 
compressible flow code has been programmed using the dual potential method. The 
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formulation (Wong and Reizes 1984) is applicable to ducts of arbitrary but constant 
cross section. The code has been used to calculate the developing flow in rectangular 
ducts of various aspect ratios for Reynolds numbers of 10 and 50. The duct geometry 
is shown in Figure 4.55. Uniform Cartesian grids were used to compute the solutions 
to be presented. 

For the three-dimensional duct geometry shown in Figure 4.55, the aspect ratio is 
b/a. White (1974) and Shah and London (1978) give a formula to compute the fully- 
developed velocity profile for constant property incompressible flow in a rectangular 
duct of constant cross-section. 

For the streamwise direction x and a rectangular section with —a < y < a and 


—b < z < b: 


u(y,z) 

Q 


16a 2 

//7T^ 

4 fea 3 
3 




g 1 
n=l,3,5,... 


192 /a\ 
~^\b) 


(_l)("-l)/2 



cosh (nnz/2a) 
cosh (nnb/2a) 


n- 


oo 1 

£ ~T tanh 

:l,3,5,... n 




Where Q is the volume flow rate, Q = ti m 4 and A = 2a x 2b. 

Note that the viscosity ft does not matter for the fully-developed velocity profile 
of The pressure gradient may be eliminated in the equation for u so that 

may be determined readily from the above equations. 

The computed flow development and fully-developed profiles agree well with 
the known results. Better agreement could be obtained using stretched grids. The 
centerline velocity development for a square duct is shown in Figure 4.56. Computed 
results using the dual potential method are indicated by symbols. The centerline 
velocity development as computed on stretched grids by Wong and Reizes (1984) 
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Table 4.2: Fully- developed centerline velocity and skin friction 



“max/ u m 

G f Re 

Aspect ratio 

Re 

dual pot. grid 

ref. 

dual pot. 

ref. 

dual pot. 

r i.oo 

10 

15 x 15 x 45 

2.096 

2.065 

14.22708 

14.07 

1.00 

50 

29 x 29 x 60 

2.096 

2.084 

14.22708 

14.23 

0.75 

50 

21 x 15 x 30 

2.077 

2.051 

14.47570 

14.38 

0.50 

50 

31 x 15 x 30 

1.992 

1.968 

15.54806 

15.40 

0.25 

50 

61 x 15 x 30 

1.774 

1.752 

18.23278 

17.90 


is shown as the Re = 10, 50 and 200 results. Experimental results of Goldstein 
and Kreid (1967) would be very near the Re = 200 results. The fully- developed 
velocity profile for Re = 50 is compared to the analytical solution in Figure 4.57 . The 
transverse velocity vectors near the exit of a square duct are shown in Figure 4.58 and 
near the exit of an aspect ratio = 0.50 rectangular duct in Figure 4.59. The velocity 
vectors are magnified 100 times to show the persistent vortical flow in the corners 
even at the fully-developed condition. A similar symmetric pattern is exhibited in all 
the rectangular ducts. The transverse flow is toward the corner along the walls and 
away from the corners along the corner bisector. Note that this flow pattern happens 
to be the opposite of the Reynolds stress driven secondary flow for a fully-developed 
turbulent flow in a rectangular duct (Demuren and Rodi 1987; Speziale 1987a). The 
flow patterns here probably result from initial disturbances that have not yet decayed. 
Stretched grids are needed for more efficient calculation of this flow. 

A summary of the fully-developed centerline velocity and wall skin friction co- 
efficient for various aspect ratio ducts is given in Table 4.2. The tolerance used in 
computing the dual potential results was 0.0001. These results were computed on a 
uniform Cartesian grid. For the square cross-section duct, the skin friction coefficient 
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Figure 4.58: Transverse velocity vectors at a plane near the exit of a square duct 

with Re = 50 
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Figure 4.59: Transverse velocity vectors at a plane near the exit of a 0.50 aspect 

ratio duct with Re = 50 


asymptotic value is computed to be C f Re = 14.23 on a 29 x 29 x 60 grid. This is the 
exact answer to two places behind the decimal point. 

The shear force was computed by integrating the product of the first order ve- 
locity derivative and the local wall area over the entire wall surface. A second order 
polynomial fit of the u velocity component is used to compute the velocity derivative 
in the shear stress equation. The local wall shear stress is computed from. 

( du dv\ 

T *v = + 

/ dw du\ 

Tiz = + 

/ dv dw\ 

T * z = ^ Val + 

The shear stresses T X y and t xz are used in this computation. The integrated wall 

shear force is then used to compute C f Re. 

The three-dimensional code has not been vectorized and uses point Gauss-Seidel 
with SOR for the vector potential Poisson equation. That particular scheme becomes 
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more and more inefficient as the number of grid points increases. The computation 
time can be reduced somewhat for steady state calculations by tightening the tol- 
erance gradually as iterations increase. For a uniform Cartesian grid, 15 x 15 x 30 
points is a good compromise of accuracy and solution speed. It has been observed 
that a solution tolerance of 0.001-0.0001 is adequate. For more accurate results on a 
uniform Cartesian grid at the Reynolds numbers reported here, the spacing as in the 
square duct case with 29 X 29 x 60 points should be used. This gives fully-developed 
flow solutions within 1% of the exact. Note that for incompressible flow the scalar 
potential can be solved once and for all for a given fixed geometry. This feature is 
used in the solution presented here. For the rectangular duct geometry, <f> — x is the 

a 

analytical solution to the Laplace equation = 0 with the boundary conditions 
discussed in Section 3.2.1. 

grid point 

on a single processor of the Cray X-MP. The complete flow field is solved by the cur- 
rent code. The symmetry of this problem was not used to speed the solution. This 
would immediately reduce the computation time by about a factor of four. The 
MFLOP rate is 12 MFLOP overall, without any enhancements. 

The three-dimensional incompressible code should be extended by adding the 
energy equation, generalized coordinates (so a curved duct case can be computed) 
and a rotating coordinate capability for centrifugal compressor modeling. A more 
efficient Poisson solver is needed and other enhancements to improve the computation 
speed. 


The computation rate for this unvectorized code is ~ i « « ■ , 

r global iteration : 
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4.4.2 Compressible flow 

The transport equations for u and B in the three-dimensional compressible 
formulation are given in Appendix A. No solutions were attempted. 
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5. DISCUSSION AND CONCLUSIONS 


A dual potential procedure has been developed and evaluated for flow calculations 
ranging from potential flow to full Navier-Stokes solutions in two dimensions. A three- 
dimensional incompressible formulation has also been presented. 

Compressibility has been handled by the dilatation variable. A dilatation trans- 
port equation is obtained from the momentum equations. Workable dilatation bound- 
ary conditions have been presented. 

From this study the following conclusions can be made: 

1. The dual potential formulation of the Navier-Stokes equations is very flexible 
in that it is easy to compute subset equations (potential, Euler). 

2. The dual potential method can simulate irrotational and rotational inflow easily. 
In fact, it is interesting to note that the test cases for which experimental results 
are available agreed best with the computed results for an irrotational inflow 
condition — boundary layer and three-dimensional duct inlet flow. Yet many 
of the computations had to be compared with the rotational inlet condition as 
this is what other computational investigators have used, namely in the channel 
heat transfer cases. Van Dyke (1970) reports that the irrotational inlet condition 
may be more realistic for experimental models with a rounded entrance. The 
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above statement and the small sampling of experimental results suggests that 
the irrotational inlet condition is an important feature to be able to simulate. 
The dual potential method handles either condition easily. 

3. Computational effort can be reduced in irrotational regions (u>= 0) and in 
incompressible regions ( B = 0). 

4. The method appears to be very accurate for incompressible calculations. 

5. It appears possible to extend the dual potential method to compressible flow, 
but the dilatation transport equation may need more work. It seems that some 
problem still persists here — perhaps self consistency of the inflow B boundary 
condition — that is most noticeable for viscous compressible problems. Also, 
the dilatation can undergo much variation in certain regions of the flow field 
making it difficult to resolve. For example, in stagnation regions B changes sign 
as shown in Figure 4.32 for the flow over a parabolic arc bump. This behavior 
of the dilatation variable may be undesirable in practical use of this method. 

6. The scalar potential requires a tight convergence tolerance as used in the test 
cases here. This may be explained by the fact that the scalar potential is the 
largest component in the velocity decomposition for the streamwise velocity. 
The convergence tolerance on the vector potential can be 1-2 orders of mag- 
nitude less restrictive than for the scalar potential. This observation may be 
problem specific since the streamwise velocity component was the dominant 
component for the test cases studied in this work. 

7. For incompressible flows, the dual potential approach can be very competitive 
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with primitive variable methods. As an example, the original INS3D code 
(Kwak et al. 1986) performs at ~ ^ grfj point on the Cray 

X-MP. The non-optimized version of the three-dimensional incompressible dual 
potential code runs at ~ gjj x grid point on lhe X ‘ MP - The onl ^ 

significant point to make about this comparison is that the dual potential code is 
written for uniform Cartesian grids while the INS3D code can handle generalized 
three-dimensional coordinates. 

8. It is possible to solve the dual potential equation set in an iterative, uncoupled 
way. 

A number of capabilities of this formulation have been demonstrated, some for 
the first time. The success seems to be limited for compressible viscous flows with a 
pressure gradient. Determination of the wall density is awkward for these cases and 
perhaps causes some problems for the dilatation transport equation. The dilatation 
transport equation seems to work fine for a known density and temperature field as 
was the case for the boundary-layer solutions. More remains to be done before a 
trouble free compressible formulation is available. 
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0. RECOMMENDATIONS FOR FUTURE WORK 


The key advantages and disadvantages of the dual potential method have been 
stated in Section 1.3. 2. 2. Speziale (1987b) describes further advantages of non- 
primitive variable methods using vorticity. He shows that these methods have simpler 
boundary conditions than primitive variable formulations for problems in a rotating 
reference frame. A review of the literature (Ozoe et al. 1985) suggests that a dual 
potential formulation is the method of choice for confined flows. The solution speed 
of the dual potential type methods for incompressible flows is well documented (Aziz 
and Heliums 1967). However, the disadvantages have precluded solutions in complex 
geometries. 

It seems best then to concentrate research efforts on the known advantages of this 
method and solve problems of practical interest. A few problems will be suggested 
that exploit the advantages of the dual potential method. Global weather forecasting 
is one such interesting practical problem. The geometry is simple, the flow is confined 
and density variations can be obtained from simple correlations. Incompressible flows 
inside curved, twisted ducts have been solved using the dual potential method (Yang 
and Camarero 1986). It should be possible to build on these solutions by allowing 
the duct to rotate. Many practical fluid flow problems could then be modeled such 
as the flow through an automotive water pump, a centrifugal compressor, a flow 
rikCObUinvi HM(aL UlMiyK IMU I rlUiVltiL) 
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meter, etc. Further interesting practical solutions could be obtained for problems that 
resemble the bump cases studied here. The “walls” in that two-dimensional problem 
can be allowed to move and the bump can have a time-dependent growth. Also, 
the flow can be treated as viscous. Since solutions have been obtained for subsonic 
compressible flows, the dual potential method may be well suited for computing flows 
with embedded compressible regions. Examples of practical flows with these features 
are the flow under an automobile and flow through a carburetor. Heat transfer effects 
can be included in all of these simulations. 

In other research areas the dual potential method may also be of importance. 
Direct simulations of channel flow and a flat plate boundary layer have been obtained 
by Rai and Moin (private communication, NASA/Ames Research Center). They 
used a non-conservative primitive variable method on a staggered grid to simulate 
incompressible turbulent flow. The dual potential method could be used to compute 
direct simulations of these flows with little modification. If the speed advantage 
of the dual potential method for incompressible flow calculations extends to direct 
simulations, then interest in the formulation will surely increase. 

It is always wise to seek out similarities in the other disciplines. The Helmholtz 
decomposition theorem was developed and used extensively by scientists working in 
electrodynamics. There may be some analogous physical problem in physics, mathe- 
matics or other field that could help in the understanding and use of the dual potential 
method. A search of the literature quickly turned up the paper by Muller (1987) on 
the topic of vector splitting applied to a physics problem. 

Since the dependent variables in the dual potential formulation are vorticity and 
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dilatation, one should use these variables to his advantage not only in computing a 
flow, but in the understanding of the physics. For example, vorticity is generated 
by such natural phenomenon as walls and shocks. A code which solves for vorticity 
directly, such as the dual potential method, may be useful in problems in which one 
wishes to track the vorticity. The dilatation variable can also be useful in understand- 
ing some flow situations. The dilatation changes sign in the vicinity of stagnation 
points and shocks. This physical behavior can help give insight into the problem 
being studied (and cause problems numerically!). 

Finally, it is important to have a code that is a good Navier-Stokes solver, but 
also is easily and efficiently used in other modes. It has been demonstrated that the 
dual potential method can be used in a potential, Euler and Navier-Stokes mode. 
Further work can improve the performance even more and address the issues of shock 
capturing, turbulence modeling, direct simulations, etc. If this method is ever to 
prove useful in applications, it must be able to run in all these modes efficiently and 
reliably. In addition it must be faster than existing methods. Unfortunately, the 
ideas behind the dual potential method are not as familiar as the primitive variable 
approaches. This unfamiliarity is a handicap since the details of application codes 
should be easily understood. To increase the familiarity of this method it must be 
tested on some problems of practical interest. 

The two- and three-dimensional dual potential codes written as a part of this 
research are still under development. A copy of a version of the codes may be obtained 
from the first author or through the Department of Mechanical Engineering, Iowa 
State University, Ames, Iowa 50010 (Attn: Professor R. H. Pletcher). 
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8. APPENDIX A: THREE-DIMENSIONAL VORTICITY AND 
DILATATION TRANSPORT EQUATIONS 


The momentum equations for a three-dimensional flow of a Newtonian Stokesian 
fluid will be presented here. They will be written in terms of vorticity (actually and 
dilatation as the dependent variables. First a point will be made about the number 
of transport equations which must be solved for the vorticity in three dimensions. 

In three dimensions it is only necessary to solve two of the three component 
vorticity transport equations because the third vorticity can always be obtained by a 
linear combination of derivatives of the other two. To see this, consider the vorticity 
definition, 

w= Vx V 

where, 

U>1 = Wy — v z 

= u z - w x 

U/% = Vx — Uy 

therefore, 

^3 2 = -“Tz ~“2y 


Pf&QtD'im PAGE BLA&K NOT FILMED 


8LANK 
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Hence, the third vorticity component can always be obtained from the other two 
components. In the above, the numbers 1,2,3 represent directions and the letters 
x,y,z are derivatives. 

The vorticity transport and the dilatation transport equations follow. They 
represent the momentum equations for a three-dimensional flow. To compute a flow 
field using the dual potential method, additional equations are required to satisfy 
mass and energy conservation. Also, the Poisson equations for the potentials must 
be solved (Equations 2.13 and 2.14): 

V 2 (f> = B 
V 2 A = -w 

Equations 8.1-8.3 display the vorticity transport equations in the dependent 
variable, . The three-dimensional dilatation transport equation is presented in 
Equation 8.4. The ideal gas law (p = pRT) has been used to replace the pressure 
gradient term in the momentum equations. Body forces are neglected. The following 
equations have not been coded. Their numerical behavior is unknown. 
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Vorticity transport: ^-component 
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Vorticity transport: y-component 
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Vorticity transport: z-component 
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Dilatation transport: 


B t + uB x + vB y + wB z - - 

— (tlx + Vy + w\ + 2 UyVx + 2 u Z W X + 2 v Z Wy ) 
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9. APPENDIX B: ALTERNATE NON-PRIMITIVE VARIABLE 

EQUATIONS 


An approach different from that presented in the main body of this report is to 
form a different combination with the x,y momentum equations. 

In two dimensions the governing equations are: 

1. continuity 

2. x momentum 

3. y momentum 

4. energy 

The cross product applied to the momentum equations yields the vorticity transport 
equation. This leaves one with still another usage of the momentum equations. In 
the main body of this report the dilatation variable, B = u x + Vy , was selected as the 
final usage. Another possible operation on the momentum equations is to compute 
(x momentum) + (y momentum). Then, choosing the dependent variable to 
be (uy+v x ), an equation for the rate of shear deformation is derived. Let T = uy+v x . 


d 

~§y 
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The two-dimensional equation for T is 

p + uT® + vTy + Ti?j = —2 p X y + j^V 2 r + p^VT ' ^ + 

1 5 2 Du Du 

+ ~ Py m ~ px m (9 ' 1) 

Then, the momentum equations provide the solution for co and I\ Recalling the 
definitions, 


u> = Vx — Uy (9-2) 

r = V® 4 -Uy (9.3) 


one obtains the following compact formulas for some velocity derivatives: 


u> + r 


r-w 

u y = — T" 


(9-4) 

(9.5) 


This adaptation to the dual potential method could be used as an inverse solution 
procedure since the wall shear stress is simply pT. 

NOTICE: At a no-slip impermeable boundary, v® = 0. There is no trouble 
getting the vorticity at the wall. Therefore, T should be simple to obtain once the 
vorticity is computed! At the wall T = -u>. Also, the potentials are not needed for 
this use of the momentum equations. 



